46 #define MIN_DEBUG_LEVEL 3
56 #define M_PI 3.14159265358979323846
60 #define GROMACS_PAIR 1
63 #ifndef GROMACS_EXCLUSIONS
64 #define GROMACS_EXCLUSIONS 1
69 #ifndef MOLECULE2_C // first object file
71 #ifdef MEM_OPT_VERSION
72 template int lookupCstPool<AtomSignature>(
const vector<AtomSignature>&,
const AtomSignature&);
73 template int lookupCstPool<ExclusionSignature>(
const vector<ExclusionSignature>&,
const ExclusionSignature&);
77 const char *segid,
int resid,
int *begin,
int *end)
const {
80 while ( elem && strcasecmp(elem->
mySegid,segid) ) elem = elem->
next;
81 if ( elem && (resid >= elem->
firstResid) && (resid <= elem->lastResid) ) {
90 const char *segid,
int resid,
int aid) {
92 if ( firstResid == -1 ) {
93 strcpy(mySegid,segid);
98 }
else if ( ! strcasecmp(mySegid,segid) ) {
99 if ( resid == lastResid ) {
100 atomIndex[lastResid - firstResid + 1] = aid + 1;
101 }
else if ( resid < lastResid ) {
104 " out of order in segment " << segid <<
105 ", lookup for additional residues in this segment disabled.\n" <<
endi;
107 next->append(segid,resid,aid);
109 for ( ; lastResid < resid; ++lastResid )
atomIndex.add(aid);
110 atomIndex[lastResid - firstResid + 1] = aid + 1;
114 next->append(segid,resid,aid);
122 const char *segid,
int resid,
const char *aname)
const {
124 if (atomNames == NULL || resLookup == NULL)
126 NAMD_die(
"Tried to find atom from name on node other than node 0");
131 if ( resLookup->lookup(segid,resid,&i,&end) )
return -1;
132 for ( ; i < end; ++i ) {
133 #ifdef MEM_OPT_VERSION
134 Index idx = atomNames[i].atomnameIdx;
137 if ( ! strcasecmp(aname,atomNames[i].atomname) )
return i;
145 const char *segid,
int resid)
const {
147 if (atomNames == NULL || resLookup == NULL)
149 NAMD_die(
"Tried to find atom from name on node other than node 0");
153 if ( resLookup->lookup(segid,resid,&i,&end) )
return 0;
159 const char *segid,
int resid,
int index)
const {
161 if (atomNames == NULL || resLookup == NULL)
163 NAMD_die(
"Tried to find atom from name on node other than node 0");
167 if ( resLookup->lookup(segid,resid,&i,&end) )
return -1;
168 if ( index >= 0 && index < ( end - i ) )
return ( index + i );
185 if (
sizeof(
int32) != 4 ) {
NAMD_bug(
"sizeof(int32) != 4"); }
187 this->params = param;
195 is_lonepairs_psf = 1;
205 lcpoParamType = NULL;
214 #ifdef MEM_OPT_VERSION
222 atomChargePool = NULL;
223 eachAtomCharge = NULL;
236 #ifndef MEM_OPT_VERSION
242 dihedralsByAtom=NULL;
243 impropersByAtom=NULL;
244 crosstermsByAtom=NULL;
246 gromacsPairByAtom=NULL;
254 #ifdef MEM_OPT_VERSION
256 exclChkSigPool = NULL;
258 eachAtomExclSig = NULL;
260 fixedAtomsSet = NULL;
261 constrainedAtomsSet = NULL;
264 fullExclusionsByAtom=NULL;
265 modExclusionsByAtom=NULL;
272 #ifdef MEM_OPT_VERSION
279 exPressureAtomFlags=NULL;
280 rigidBondLengths=NULL;
295 consTorqueIndexes=NULL;
296 consTorqueParams=NULL;
297 consForceIndexes=NULL;
301 alch_unpert_bonds = NULL;
302 alch_unpert_angles = NULL;
303 alch_unpert_dihedrals = NULL;
313 #ifndef MEM_OPT_VERSION
351 numExPressureAtoms=0;
353 numFixedRigidBonds=0;
354 numMultipleDihedrals=0;
355 numMultipleImpropers=0;
364 numCalcFullExclusions=0;
372 num_alch_unpert_Bonds = 0;
373 num_alch_unpert_Angles = 0;
374 num_alch_unpert_Dihedrals = 0;
431 initialize(simParams,param);
444 initialize(simParams,param);
446 #ifdef MEM_OPT_VERSION
448 read_mol_signatures(filename, param, cfgList);
450 read_psf_file(filename, param);
453 assignLCPOTypes( 0 );
466 #ifdef MEM_OPT_VERSION
467 NAMD_die(
"Sorry, plugin IO is not supported in the memory optimized version.");
469 initialize(simParams, param);
471 int optflags = MOLFILE_BADOPTIONS;
472 molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*
sizeof(molfile_atom_t));
473 memset(atomarray, 0, natoms*
sizeof(molfile_atom_t));
476 int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
477 if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
479 NAMD_die(
"ERROR: plugin failed reading structure data");
481 if(optflags == MOLFILE_BADOPTIONS) {
483 NAMD_die(
"ERROR: plugin didn't initialize optional data flags");
485 if(optflags & MOLFILE_OCCUPANCY) {
486 setOccupancyData(atomarray);
488 if(optflags & MOLFILE_BFACTOR) {
489 setBFactorData(atomarray);
492 plgLoadAtomBasics(atomarray);
499 int *bondtype, nbondtypes;
501 if(pIOHdl->read_bonds!=NULL) {
502 if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder,
503 &bondtype, &nbondtypes, &bondtypename)){
504 NAMD_die(
"ERROR: failed reading bond information.");
509 plgLoadBonds(from,to);
513 int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
514 int ctermcols, ctermrows;
515 int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes;
516 int *impropertypes, numimpropertypes;
517 char **angletypenames, **dihedraltypenames, **impropertypenames;
519 plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
520 if(pIOHdl->read_angles!=NULL) {
521 if(pIOHdl->read_angles(pIOFileHdl,
522 &numAngles, &plgAngles,
523 &angletypes, &numangletypes, &angletypenames,
524 &numDihedrals, &plgDihedrals,
525 &dihedraltypes, &numdihedraltypes, &dihedraltypenames,
526 &numImpropers, &plgImpropers,
527 &impropertypes, &numimpropertypes, &impropertypenames,
528 &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) {
529 NAMD_die(
"ERROR: failed reading angle information.");
533 if(numAngles!=0) plgLoadAngles(plgAngles);
534 if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
535 if(numImpropers!=0) plgLoadImpropers(plgImpropers);
536 if(numCrossterms!=0) plgLoadCrossterms(plgCterms);
538 numRealBonds = numBonds;
542 assignLCPOTypes( 2 );
565 if (atomNames != NULL)
572 if (resLookup != NULL)
576 if (drudeConsts != NULL)
delete [] drudeConsts;
577 if (lphosts != NULL)
delete [] lphosts;
578 if (anisos != NULL)
delete [] anisos;
579 if (tholes != NULL)
delete [] tholes;
580 if (lphostIndexes != NULL)
delete [] lphostIndexes;
584 if (lcpoParamType != NULL)
delete [] lcpoParamType;
586 #ifdef MEM_OPT_VERSION
587 if(eachAtomSig)
delete [] eachAtomSig;
596 if (dihedrals != NULL)
599 if (impropers != NULL)
602 if (crossterms != NULL)
603 delete [] crossterms;
612 if (acceptors != NULL)
615 #ifdef MEM_OPT_VERSION
616 if(exclSigPool)
delete [] exclSigPool;
617 if(exclChkSigPool)
delete [] exclChkSigPool;
618 if(eachAtomExclSig)
delete [] eachAtomExclSig;
619 if(fixedAtomsSet)
delete fixedAtomsSet;
620 if(constrainedAtomsSet)
delete constrainedAtomsSet;
622 if (bondsByAtom != NULL)
623 delete [] bondsByAtom;
625 if (anglesByAtom != NULL)
626 delete [] anglesByAtom;
628 if (dihedralsByAtom != NULL)
629 delete [] dihedralsByAtom;
631 if (impropersByAtom != NULL)
632 delete [] impropersByAtom;
634 if (crosstermsByAtom != NULL)
635 delete [] crosstermsByAtom;
640 if (fullExclusionsByAtom != NULL)
641 delete [] fullExclusionsByAtom;
643 if (modExclusionsByAtom != NULL)
644 delete [] modExclusionsByAtom;
646 if (all_exclusions != NULL)
647 delete [] all_exclusions;
650 if (gromacsPairByAtom != NULL)
651 delete [] gromacsPairByAtom;
655 if (tholesByAtom != NULL)
656 delete [] tholesByAtom;
657 if (anisosByAtom != NULL)
658 delete [] anisosByAtom;
663 if (lcpoParamType != NULL)
664 delete [] lcpoParamType;
666 if (fixedAtomFlags != NULL)
667 delete [] fixedAtomFlags;
669 if (stirIndexes != NULL)
670 delete [] stirIndexes;
673 #ifdef MEM_OPT_VERSION
674 if(clusterSigs != NULL){
675 delete [] clusterSigs;
681 if (clusterSize != NULL)
682 delete [] clusterSize;
684 if (exPressureAtomFlags != NULL)
685 delete [] exPressureAtomFlags;
687 if (rigidBondLengths != NULL)
688 delete [] rigidBondLengths;
691 if (fepAtomFlags != NULL)
692 delete [] fepAtomFlags;
693 if (alch_unpert_bonds != NULL)
694 delete [] alch_unpert_bonds;
695 if (alch_unpert_angles != NULL)
696 delete [] alch_unpert_angles;
697 if (alch_unpert_dihedrals != NULL)
698 delete [] alch_unpert_dihedrals;
702 if (ssAtomFlags != NULL)
703 delete [] ssAtomFlags;
704 if (ss_vdw_type != NULL)
705 delete [] ss_vdw_type;
706 if (ss_index != NULL)
710 if (qmAtomGroup != NULL)
711 delete [] qmAtomGroup;
713 if (qmAtmIndx != NULL)
716 if (qmAtmChrg != NULL)
720 if (qmGrpNumBonds != NULL)
721 delete [] qmGrpNumBonds;
723 if (qmGrpSizes != NULL)
724 delete [] qmGrpSizes;
726 if (qmDummyBondVal != NULL)
727 delete [] qmDummyBondVal;
729 if (qmMMNumTargs != NULL)
730 delete [] qmMMNumTargs;
735 if (qmGrpChrg != NULL)
738 if (qmGrpMult != NULL)
741 if (qmMeMMindx != NULL)
742 delete [] qmMeMMindx;
744 if (qmMeQMGrp != NULL)
747 if (qmLSSSize != NULL)
750 if (qmLSSIdxs != NULL)
753 if (qmLSSMass != NULL)
756 if (qmLSSRefSize != NULL)
757 delete [] qmLSSRefSize;
759 if (qmLSSRefIDs != NULL)
760 delete [] qmLSSRefIDs;
762 if (qmLSSRefMass != NULL)
763 delete [] qmLSSRefMass;
765 if (qmMMBond != NULL) {
766 for(
int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
767 if (qmMMBond[grpIndx] != NULL)
768 delete [] qmMMBond[grpIndx];
773 if (qmGrpBonds != NULL) {
774 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
775 if (qmGrpBonds[grpIndx] != NULL)
776 delete [] qmGrpBonds[grpIndx];
778 delete [] qmGrpBonds;
781 if (qmMMBondedIndx != NULL) {
782 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
783 if (qmMMBondedIndx[grpIndx] != NULL)
784 delete [] qmMMBondedIndx[grpIndx];
786 delete [] qmMMBondedIndx;
789 if (qmMMChargeTarget != NULL) {
790 for(
int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
791 if (qmMMChargeTarget[grpIndx] != NULL)
792 delete [] qmMMChargeTarget[grpIndx];
794 delete [] qmMMChargeTarget;
797 if (qmElementArray != NULL){
798 for(
int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
799 if (qmElementArray[atmInd] != NULL)
800 delete [] qmElementArray[atmInd];
802 delete [] qmElementArray;
805 if (qmDummyElement != NULL){
806 for(
int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
807 if (qmDummyElement[atmInd] != NULL)
808 delete [] qmDummyElement[atmInd];
810 delete [] qmDummyElement;
813 if (qmCustPCSizes != NULL){
814 delete [] qmCustPCSizes;
817 if (qmCustomPCIdxs != NULL){
818 delete [] qmCustomPCIdxs;
821 if (cSMDindex != NULL) {
822 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
823 if (cSMDindex[grpIndx] != NULL)
824 delete [] cSMDindex[grpIndx];
829 if (cSMDpairs != NULL) {
830 for(
int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) {
831 if (cSMDpairs[instIndx] != NULL)
832 delete [] cSMDpairs[instIndx];
837 if (cSMDindxLen != NULL)
838 delete [] cSMDindxLen;
841 if (cSMDVels != NULL)
843 if (cSMDcoffs != NULL)
846 #ifndef MEM_OPT_VERSION
853 #ifndef MEM_OPT_VERSION
874 void Molecule::read_psf_file(
char *fname,
Parameters *params)
884 if ( (psf_file =
Fopen(fname,
"r")) == NULL)
886 sprintf(err_msg,
"UNABLE TO OPEN .psf FILE %s", fname);
902 sprintf(err_msg,
"EMPTY .psf FILE %s", fname);
910 sprintf(err_msg,
"UNABLE TO FIND \"PSF\" STRING IN PSF FILE %s",
919 iout <<
iWARN <<
"Reading PSF supporting DRUDE without "
920 "enabling the Drude model in the simulation config file\n" <<
endi;
939 sprintf(err_msg,
"MISSING EVERYTHING BUT PSF FROM %s", fname);
947 sprintf(err_msg,
"CAN NOT FIND \"NTITLE\" STRING IN PSF FILE %s",
952 sscanf(buffer,
"%d", &NumTitle);
968 sprintf(err_msg,
"FOUND EOF INSTEAD OF NATOM IN PSF FILE %s",
981 sprintf(err_msg,
"DIDN'T FIND \"NATOM\" IN PSF FILE %s",
987 sscanf(buffer,
"%d", &numAtoms);
989 read_atoms(psf_file, params);
1002 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NBONDS IN PSF");
1008 NAMD_die(
"DID NOT FIND NBOND AFTER ATOM LIST IN PSF");
1012 sscanf(buffer,
"%d", &numBonds);
1015 read_bonds(psf_file);
1028 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NTHETA IN PSF");
1034 NAMD_die(
"DID NOT FIND NTHETA AFTER BOND LIST IN PSF");
1038 sscanf(buffer,
"%d", &numAngles);
1041 read_angles(psf_file, params);
1054 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NPHI IN PSF");
1060 NAMD_die(
"DID NOT FIND NPHI AFTER ANGLE LIST IN PSF");
1064 sscanf(buffer,
"%d", &numDihedrals);
1067 read_dihedrals(psf_file, params);
1080 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NIMPHI IN PSF");
1086 NAMD_die(
"DID NOT FIND NIMPHI AFTER ATOM LIST IN PSF");
1090 sscanf(buffer,
"%d", &numImpropers);
1093 read_impropers(psf_file, params);
1106 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NDON IN PSF");
1112 NAMD_die(
"DID NOT FIND NDON AFTER ATOM LIST IN PSF");
1116 sscanf(buffer,
"%d", &numDonors);
1119 read_donors(psf_file);
1132 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NACC IN PSF");
1138 NAMD_die(
"DID NOT FIND NACC AFTER ATOM LIST IN PSF");
1142 sscanf(buffer,
"%d", &numAcceptors);
1145 read_acceptors(psf_file);
1154 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NNB IN PSF FILE");
1159 sscanf(buffer,
"%d", &numExclusions);
1162 read_exclusions(psf_file);
1181 int is_found_numlp = 0;
1182 int is_found_numaniso = 0;
1183 int is_found_ncrterm = 0;
1184 while ( ! is_found ) {
1190 if (ret_code != 0) {
1193 else if ( (is_found_numlp =
NAMD_find_word(buffer,
"NUMLP")) > 0) {
1194 is_found = is_found_numlp;
1196 else if ( (is_found_numaniso =
NAMD_find_word(buffer,
"NUMANISO")) > 0) {
1197 is_found = is_found_numaniso;
1199 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1200 is_found = is_found_ncrterm;
1204 if (is_found_numlp) {
1207 sscanf(buffer,
"%d", &numLphosts);
1210 if (numLphosts == 0) {
1215 is_lonepairs_psf = 0;
1221 NAMD_die(
"FOUND LONE PAIR HOSTS IN PSF WITH \"LONEPAIRS\" DISABLED IN CONFIG FILE");
1225 if (numBonds) process_bonds(params);
1227 if (is_found_numlp) {
1229 if (numLphosts > 0) read_lphosts(psf_file);
1233 is_found_numaniso = 0;
1234 is_found_ncrterm = 0;
1235 while ( ! is_found ) {
1241 if (ret_code != 0) {
1244 else if ( (is_found_numaniso =
NAMD_find_word(buffer,
"NUMANISO")) > 0) {
1245 is_found = is_found_numaniso;
1247 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1248 is_found = is_found_ncrterm;
1253 if (is_found_numaniso && is_drude_psf) {
1256 sscanf(buffer,
"%d", &numAnisos);
1257 if (numAnisos > 0) read_anisos(psf_file);
1261 is_found_ncrterm = 0;
1262 while ( ! is_found ) {
1268 if (ret_code != 0) {
1271 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1272 is_found = is_found_ncrterm;
1276 else if (is_drude_psf ) {
1277 NAMD_die(
"DID NOT FIND REQUIRED NUMANISO IN DRUDE PSF FILE");
1279 else if (is_found_numaniso ) {
1280 NAMD_die(
"FOUND NUMANISO IN PSF FILE MISSING DRUDE DESIGNATION");
1283 if (is_found_ncrterm) {
1286 sscanf(buffer,
"%d", &numCrossterms);
1287 if (numCrossterms > 0) read_crossterms(psf_file, params);
1295 if (is_drude_psf && numDrudeAtoms) {
1297 Bond *newbonds =
new Bond[numBonds+numDrudeAtoms];
1298 memcpy(newbonds, bonds, numBonds*
sizeof(
Bond));
1301 int origNumBonds = numBonds;
1302 for (i=0; i < numAtoms; i++) {
1303 if (!is_drude(i))
continue;
1304 Bond *b = &(bonds[numBonds++]);
1308 atomNames[i-1].atomtype, atomNames[i].atomtype, b
1311 if (numBonds-origNumBonds != numDrudeAtoms) {
1312 NAMD_die(
"must have same number of Drude particles and parents");
1317 numRealBonds = numBonds;
1318 build_atom_status();
1339 void Molecule::read_atoms(FILE *fd,
Parameters *params)
1344 int last_atom_number=0;
1346 char segment_name[11];
1347 char residue_number[11];
1348 char residue_name[11];
1368 hydrogenGroup.resize(0);
1370 if (
atoms == NULL || atomNames == NULL )
1372 NAMD_die(
"memory allocation failed in Molecule::read_atoms");
1378 while (atom_number < numAtoms)
1391 read_count=sscanf(buffer,
"%d %s %s %s %s %s %f %f",
1392 &atom_number, segment_name, residue_number,
1393 residue_name, atom_name, atom_type, &charge, &mass);
1394 if (mass <= 0.05) ++numZeroMassAtoms;
1397 if (read_count != 8)
1401 sprintf(err_msg,
"BAD ATOM LINE FORMAT IN PSF FILE IN ATOM LINE %d\nLINE=%s",
1402 last_atom_number+1, buffer);
1415 read_count=sscanf(buffer,
1418 "%*d %*s %*s %*s %*s %*s %*f %*f %*d %f %f", &alpha, &thole);
1419 if (read_count != 2)
1423 sprintf(err_msg,
"BAD ATOM LINE FORMAT IN PSF FILE "
1424 "IN ATOM LINE %d\nLINE=%s", last_atom_number+1, buffer);
1427 drudeConsts[atom_number-1].alpha = alpha;
1428 drudeConsts[atom_number-1].thole = thole;
1429 if (fabs(alpha) >= 1e-6) ++numDrudeAtoms;
1435 if ( sscanf(atom_type,
"%d", &atom_type_num) > 0 )
1437 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.");
1441 if (atom_number != last_atom_number+1)
1445 sprintf(err_msg,
"ATOM NUMBERS OUT OF ORDER AT ATOM #%d OF PSF FILE",
1446 last_atom_number+1);
1455 int reslength = strlen(residue_name)+1;
1456 int namelength = strlen(atom_name)+1;
1457 int typelength = strlen(atom_type)+1;
1459 atomNames[atom_number-1].resname = nameArena->getNewArray(reslength);
1460 atomNames[atom_number-1].atomname = nameArena->getNewArray(namelength);
1461 atomNames[atom_number-1].atomtype = nameArena->getNewArray(typelength);
1463 if (atomNames[atom_number-1].resname == NULL)
1465 NAMD_die(
"memory allocation failed in Molecule::read_atoms");
1469 strcpy(atomNames[atom_number-1].resname, residue_name);
1470 strcpy(atomNames[atom_number-1].atomname, atom_name);
1471 strcpy(atomNames[atom_number-1].atomtype, atom_type);
1472 atoms[atom_number-1].mass = mass;
1477 if ( tmpResLookup ) tmpResLookup =
1478 tmpResLookup->
append(segment_name, atoi(residue_number), atom_number-1);
1482 memcpy(one->
segname, segment_name, strlen(segment_name)+1);
1483 one->
resid = atoi(residue_number);
1488 }
else if (
atoms[atom_number-1].mass <= 0.05) {
1490 }
else if (
atoms[atom_number-1].mass < 1.0) {
1492 }
else if (
atoms[atom_number-1].mass <=3.5) {
1494 }
else if ((atomNames[atom_number-1].atomname[0] ==
'O') &&
1495 (
atoms[atom_number-1].mass >= 14.0) &&
1496 (
atoms[atom_number-1].mass <= 18.0)) {
1502 &(
atoms[atom_number-1]));
1521 void Molecule::read_bonds(FILE *fd)
1529 bonds=
new Bond[numBonds];
1533 NAMD_die(
"memory allocations failed in Molecule::read_bonds");
1537 while (num_read < numBonds)
1549 if (atom_nums[j] >= numAtoms)
1553 sprintf(err_msg,
"BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1559 Bond *b = &(bonds[num_read]);
1560 b->
atom1=atom_nums[0];
1561 b->
atom2=atom_nums[1];
1570 void Molecule::process_bonds(
Parameters *params) {
1573 int origNumBonds = numBonds;
1575 int numZeroFrcBonds = 0;
1577 int numDrudeBonds = 0;
1579 for (
int old_read=0; old_read < origNumBonds; ++old_read)
1582 Bond *b = &(bonds[num_read]);
1583 *b = bonds[old_read];
1584 atom_nums[0] = b->
atom1;
1585 atom_nums[1] = b->
atom2;
1590 atomNames[atom_nums[0]].atomtype,
1591 atomNames[atom_nums[1]].atomtype,
1600 numZeroFrcBonds += (k == 0.);
1601 numLPBonds += is_lp_bond;
1602 numDrudeBonds += is_drude_bond;
1612 if ( is_lonepairs_psf ) {
1613 if (k == 0. || is_lp_bond) --numBonds;
1616 numLPBonds -= is_lp_bond;
1617 if (k == 0. && !is_lp_bond) --numBonds;
1626 if (k == 0. || is_lp_bond || is_drude_bond) --numBonds;
1632 if ( num_read != numBonds ) {
1633 NAMD_bug(
"num_read != numBonds in Molecule::process_bonds()");
1637 if ( numBonds != origNumBonds ) {
1638 if (numZeroFrcBonds) {
1639 iout <<
iWARN <<
"Ignored " << numZeroFrcBonds <<
1640 " bonds with zero force constants.\n" <<
1641 iWARN <<
"Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
1645 iout <<
iWARN <<
"Ignored " << numLPBonds <<
1646 " bonds with lone pairs.\n" <<
1647 iWARN <<
"Will infer lonepair bonds from LPhost entries.\n" <<
endi;
1649 if (numDrudeBonds) {
1650 iout <<
iWARN <<
"Ignored " << numDrudeBonds <<
1651 " bonds with Drude particles.\n" <<
1652 iWARN <<
"Will use polarizability to assign Drude bonds.\n" <<
endi;
1665 alch_unpert_bonds=
new Bond[num_alch_unpert_Bonds];
1667 if (alch_unpert_bonds == NULL) {
1668 NAMD_die(
"memory allocations failed in Molecule::read_alch_unpert_bonds");
1671 while (num_read < num_alch_unpert_Bonds) {
1672 for (j=0; j<2; j++) {
1675 if (atom_nums[j] >= numAtoms) {
1678 sprintf(err_msg,
"BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN ALCH PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1683 Bond *b = &(alch_unpert_bonds[num_read]);
1684 b->
atom1=atom_nums[0];
1685 b->
atom2=atom_nums[1];
1709 void Molecule::read_angles(FILE *fd,
Parameters *params)
1715 int origNumAngles = numAngles;
1717 angles=
new Angle[numAngles];
1721 NAMD_die(
"memory allocation failed in Molecule::read_angles");
1725 while (num_read < numAngles)
1738 if (atom_nums[j] >= numAtoms)
1742 sprintf(err_msg,
"ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1748 angles[num_read].atom1=atom_nums[0];
1749 angles[num_read].atom2=atom_nums[1];
1750 angles[num_read].atom3=atom_nums[2];
1755 atomNames[atom_nums[0]].atomtype,
1756 atomNames[atom_nums[1]].atomtype,
1757 atomNames[atom_nums[2]].atomtype,
1758 &(angles[num_read]), simParams->
alchOn ? -1 : 0);
1759 if ( angles[num_read].angle_type == -1 ) {
1760 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n"
1765 Real k, t0, k_ub, r_ub;
1766 if ( angles[num_read].angle_type == -1 ) { k = -1.; k_ub = -1.; }
else
1768 if ( k == 0. && k_ub == 0. ) --numAngles;
1773 if ( numAngles != origNumAngles ) {
1774 iout <<
iWARN <<
"Ignored " << origNumAngles - numAngles <<
1775 " angles with zero force constants.\n" <<
endi;
1788 alch_unpert_angles=
new Angle[num_alch_unpert_Angles];
1790 if (alch_unpert_angles == NULL) {
1791 NAMD_die(
"memory allocation failed in Molecule::read_alch_unpert_angles");
1794 while (num_read < num_alch_unpert_Angles) {
1795 for (j=0; j<3; j++) {
1798 if (atom_nums[j] >= numAtoms) {
1800 sprintf(err_msg,
"ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN ALCH UNPERT PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1805 alch_unpert_angles[num_read].atom1=atom_nums[0];
1806 alch_unpert_angles[num_read].atom2=atom_nums[1];
1807 alch_unpert_angles[num_read].atom3=atom_nums[2];
1830 void Molecule::read_dihedrals(FILE *fd,
Parameters *params)
1833 int last_atom_nums[4];
1837 Bool duplicate_bond;
1842 last_atom_nums[j] = -1;
1845 dihedrals =
new Dihedral[numDihedrals];
1847 if (dihedrals == NULL)
1849 NAMD_die(
"memory allocation failed in Molecule::read_dihedrals");
1853 while (num_read < numDihedrals)
1855 duplicate_bond =
TRUE;
1867 if (atom_nums[j] >= numAtoms)
1871 sprintf(err_msg,
"DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1876 if (atom_nums[j] != last_atom_nums[j])
1878 duplicate_bond =
FALSE;
1881 last_atom_nums[j] = atom_nums[j];
1891 if (multiplicity == 2)
1893 numMultipleDihedrals++;
1903 dihedrals[num_unique-1].atom1=atom_nums[0];
1904 dihedrals[num_unique-1].atom2=atom_nums[1];
1905 dihedrals[num_unique-1].atom3=atom_nums[2];
1906 dihedrals[num_unique-1].atom4=atom_nums[3];
1910 atomNames[atom_nums[0]].atomtype,
1911 atomNames[atom_nums[1]].atomtype,
1912 atomNames[atom_nums[2]].atomtype,
1913 atomNames[atom_nums[3]].atomtype,
1914 &(dihedrals[num_unique-1]),
1915 multiplicity, simParams->
alchOn ? -1 : 0);
1916 if ( dihedrals[num_unique-1].dihedral_type == -1 ) {
1917 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n"
1924 numDihedrals = num_unique;
1936 alch_unpert_dihedrals =
new Dihedral[num_alch_unpert_Dihedrals];
1938 if (alch_unpert_dihedrals == NULL) {
1939 NAMD_die(
"memory allocation failed in Molecule::read_alch_unpert_dihedrals");
1942 while (num_read < num_alch_unpert_Dihedrals) {
1943 for (j=0; j<4; j++) {
1946 if (atom_nums[j] >= numAtoms) {
1949 sprintf(err_msg,
"DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN ALCH UNPERT PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1954 alch_unpert_dihedrals[num_read].atom1=atom_nums[0];
1955 alch_unpert_dihedrals[num_read].atom2=atom_nums[1];
1956 alch_unpert_dihedrals[num_read].atom3=atom_nums[2];
1957 alch_unpert_dihedrals[num_read].atom4=atom_nums[3];
1980 void Molecule::read_impropers(FILE *fd,
Parameters *params)
1983 int last_atom_nums[4];
1987 Bool duplicate_bond;
1993 last_atom_nums[j] = -1;
1996 impropers=
new Improper[numImpropers];
1998 if (impropers == NULL)
2000 NAMD_die(
"memory allocation failed in Molecule::read_impropers");
2004 while (num_read < numImpropers)
2006 duplicate_bond =
TRUE;
2018 if (atom_nums[j] >= numAtoms)
2022 sprintf(err_msg,
"IMPROPERS INDEX %d GREATER THAN NATOM %d IN IMPROPERS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
2026 if (atom_nums[j] != last_atom_nums[j])
2028 duplicate_bond =
FALSE;
2031 last_atom_nums[j] = atom_nums[j];
2042 if (multiplicity == 2)
2045 numMultipleImpropers++;
2056 impropers[num_unique-1].atom1=atom_nums[0];
2057 impropers[num_unique-1].atom2=atom_nums[1];
2058 impropers[num_unique-1].atom3=atom_nums[2];
2059 impropers[num_unique-1].atom4=atom_nums[3];
2063 atomNames[atom_nums[0]].atomtype,
2064 atomNames[atom_nums[1]].atomtype,
2065 atomNames[atom_nums[2]].atomtype,
2066 atomNames[atom_nums[3]].atomtype,
2067 &(impropers[num_unique-1]),
2076 numImpropers = num_unique;
2096 void Molecule::read_crossterms(FILE *fd,
Parameters *params)
2100 int last_atom_nums[8];
2103 Bool duplicate_bond;
2108 last_atom_nums[j] = -1;
2111 crossterms=
new Crossterm[numCrossterms];
2113 if (crossterms == NULL)
2115 NAMD_die(
"memory allocation failed in Molecule::read_crossterms");
2119 while (num_read < numCrossterms)
2121 duplicate_bond =
TRUE;
2133 if (atom_nums[j] >= numAtoms)
2137 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);
2141 if (atom_nums[j] != last_atom_nums[j])
2143 duplicate_bond =
FALSE;
2146 last_atom_nums[j] = atom_nums[j];
2152 iout <<
iWARN <<
"Duplicate cross-term detected.\n" <<
endi;
2156 crossterms[num_read].atom1=atom_nums[0];
2157 crossterms[num_read].atom2=atom_nums[1];
2158 crossterms[num_read].atom3=atom_nums[2];
2159 crossterms[num_read].atom4=atom_nums[3];
2160 crossterms[num_read].atom5=atom_nums[4];
2161 crossterms[num_read].atom6=atom_nums[5];
2162 crossterms[num_read].atom7=atom_nums[6];
2163 crossterms[num_read].atom8=atom_nums[7];
2167 atomNames[atom_nums[0]].atomtype,
2168 atomNames[atom_nums[1]].atomtype,
2169 atomNames[atom_nums[2]].atomtype,
2170 atomNames[atom_nums[3]].atomtype,
2171 atomNames[atom_nums[4]].atomtype,
2172 atomNames[atom_nums[5]].atomtype,
2173 atomNames[atom_nums[6]].atomtype,
2174 atomNames[atom_nums[7]].atomtype,
2175 &(crossterms[num_read]));
2177 if(!duplicate_bond) num_read++;
2180 numCrossterms = num_read;
2203 void Molecule::read_donors(FILE *fd)
2211 donors=
new Bond[numDonors];
2215 NAMD_die(
"memory allocations failed in Molecule::read_donors");
2219 while (num_read < numDonors)
2231 if (d[j] >= numAtoms)
2236 "DONOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE",
2237 d[j]+1, numAtoms, num_read+1);
2247 Bond *b = &(donors[num_read]);
2277 void Molecule::read_acceptors(FILE *fd)
2285 acceptors=
new Bond[numAcceptors];
2287 if (acceptors == NULL)
2289 NAMD_die(
"memory allocations failed in Molecule::read_acceptors");
2293 while (num_read < numAcceptors)
2305 if (d[j] >= numAtoms)
2309 sprintf(err_msg,
"ACCEPTOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE", d[j]+1, numAtoms, num_read+1);
2319 Bond *b = &(acceptors[num_read]);
2358 void Molecule::read_exclusions(FILE *fd)
2361 int *exclusion_atoms;
2362 register int num_read=0;
2365 register int insert_index=0;
2370 exclusion_atoms =
new int[numExclusions];
2372 if ( (
exclusions == NULL) || (exclusion_atoms == NULL) )
2374 NAMD_die(
"memory allocation failed in Molecule::read_exclusions");
2378 for (num_read=0; num_read<numExclusions; num_read++)
2386 if (exclusion_atoms[num_read] >= numAtoms)
2390 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN PSF FILE", exclusion_atoms[num_read]+1, numAtoms, num_read+1);
2399 for (num_read=0; num_read<numAtoms; num_read++)
2405 if (current_index>numExclusions)
2409 sprintf(err_msg,
"EXCLUSION INDEX %d LARGER THAN NUMBER OF EXLCUSIONS %d IN PSF FILE, EXCLUSION #%d\n",
2410 current_index+1, numExclusions, num_read);
2417 if (current_index != last_index)
2423 for (insert_index=last_index;
2424 insert_index<current_index; insert_index++)
2431 int a2 = exclusion_atoms[insert_index];
2435 }
else if ( a2 < a1 ) {
2440 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN PSF FILE\n", a1+1);
2445 last_index=current_index;
2450 delete [] exclusion_atoms;
2466 void Molecule::read_exclusions(
int* atom_i,
int* atom_j,
int num_exclusion)
2471 int loop_counter = 0;
2477 NAMD_die(
"memory allocation failed in Molecule::read_exclusions");
2481 for (loop_counter = 0; loop_counter < num_exclusion; loop_counter++) {
2483 if ( (atom_i == NULL) || (atom_j == NULL) ) {
2484 NAMD_die(
"null pointer expection in Molecule::read_exclusions");
2487 a = atom_i[loop_counter];
2488 b = atom_j[loop_counter];
2500 iout <<
iINFO <<
"ADDED " << num_exclusion <<
" EXPLICIT EXCLUSIONS: THIS VALUE WILL *NOT* BE ADDED TO THE STRUCTURE SUMMARY\n" <<
endi;
2524 void Molecule::read_lphosts(FILE *fd)
2536 Real value1, value2, value3;
2539 lphosts =
new Lphost[numLphosts];
2540 if (lphosts == NULL) {
2541 NAMD_die(
"memory allocation failed in Molecule::read_lphosts");
2543 for (i = 0; i < numLphosts; i++) {
2546 read_count=sscanf(buffer,
"%d %d %6s %f %f %f",
2547 &numhosts, &index, weight, &value1, &value2, &value3);
2552 if (read_count != 6 || index != next_index ||
2553 strcmp(weight,
"F") != 0 || numhosts < 2 || 3 < numhosts) {
2555 sprintf(err_msg,
"BAD FORMAT FOR LPHOST LINE %d IN PSF FILE LINE\n"
2556 "LINE=%s\n", i+1, buffer);
2559 lphosts[i].numhosts = numhosts;
2560 next_index += numhosts + 1;
2561 if (numhosts == 2) {
2562 lphosts[i].distance = value1;
2563 lphosts[i].angle = value2;
2564 lphosts[i].dihedral = 0.0;
2567 lphosts[i].distance = value1;
2568 lphosts[i].angle = value2 * (
M_PI/180);
2569 lphosts[i].dihedral = value3 * (
M_PI/180);
2573 Bond *newbonds =
new Bond[numBonds+numLphosts];
2574 memcpy(newbonds, bonds, numBonds*
sizeof(
Bond));
2577 for (i = 0; i < numLphosts; i++) {
2583 lphosts[i].atom4 = ( lphosts[i].numhosts == 3 ?
2586 Bond *b = &(bonds[numBonds++]);
2587 b->
atom1 = lphosts[i].atom1;
2588 b->
atom2 = lphosts[i].atom2;
2603 void Molecule::read_anisos(FILE *fd)
2606 int numhosts, index, i, read_count;
2609 anisos =
new Aniso[numAnisos];
2612 NAMD_die(
"memory allocation failed in Molecule::read_anisos");
2614 for (i = 0; i < numAnisos; i++)
2618 read_count=sscanf(buffer,
"%f %f %f", &k11, &k22, &k33);
2619 if (read_count != 3)
2622 sprintf(err_msg,
"BAD FORMAT FOR ANISO LINE %d IN PSF FILE LINE\n"
2623 "LINE=%s\n", i+1, buffer);
2626 anisos[i].k11 = k11;
2627 anisos[i].k22 = k22;
2628 anisos[i].k33 = k33;
2630 for (i = 0; i < numAnisos; i++) {
2643 if (atomType[0] ==
'H' || atomType[0] ==
'h') {
2647 }
else if (atomType[0] ==
'C' || atomType[0] ==
'c') {
2650 strcmp(atomType,
"CT" )==0 )
2658 else if (numBonds == 2)
2660 else if (numBonds == 3)
2662 else if (numBonds == 4)
2670 else if (numBonds == 3)
2677 }
else if (atomType[0] ==
'N' || atomType[0] ==
'n') {
2678 if ( strcmp(atomType,
"N3" ) == 0 ) {
2681 else if (numBonds == 2)
2683 else if (numBonds == 3)
2691 else if (numBonds == 2)
2693 else if (numBonds == 3)
2700 }
else if (atomType[0] ==
'O' || atomType[0] ==
'o') {
2702 if ( strcmp(atomType,
"O" )==0) {
2704 }
else if (strcmp(atomType,
"O2" )==0) {
2709 else if (numBonds == 2)
2716 }
else if (atomType[0] ==
'S' || atomType[0] ==
's') {
2717 if ( strcmp(atomType,
"SH" )==0) {
2724 }
else if (atomType[0] ==
'P' || atomType[0] ==
'p') {
2727 else if (numBonds == 4)
2731 }
else if (atomType[0] ==
'Z') {
2733 }
else if ( strcmp(atomType,
"MG" )==0) {
2744 if (atomType[0] ==
'H') {
2748 }
else if (atomType[0] ==
'C') {
2750 atomType[1] ==
'T' ||
2751 strcmp(atomType,
"CP1" )==0 ||
2752 strcmp(atomType,
"CP2" )==0 ||
2753 strcmp(atomType,
"CP3" )==0 ||
2754 strcmp(atomType,
"CS" )==0 ) {
2757 else if (numBonds == 2)
2759 else if (numBonds == 3)
2761 else if (numBonds == 4)
2767 strcmp(atomType,
"C" )==0 ||
2768 strcmp(atomType,
"CA" )==0 ||
2769 strcmp(atomType,
"CC" )==0 ||
2770 strcmp(atomType,
"CD" )==0 ||
2771 strcmp(atomType,
"CN" )==0 ||
2772 strcmp(atomType,
"CY" )==0 ||
2773 strcmp(atomType,
"C3" )==0 ||
2774 strcmp(atomType,
"CE1" )==0 ||
2775 strcmp(atomType,
"CE2" )==0 ||
2776 strcmp(atomType,
"CST" )==0 ||
2777 strcmp(atomType,
"CAP" )==0 ||
2778 strcmp(atomType,
"COA" )==0 ||
2779 strcmp(atomType,
"CPT" )==0 ||
2780 strcmp(atomType,
"CPH1")==0 ||
2781 strcmp(atomType,
"CPH2")==0
2785 else if (numBonds == 3)
2794 }
else if (atomType[0] ==
'N') {
2799 strcmp(atomType,
"NH3" )==0 ||
2802 strcmp(atomType,
"NP" )==0
2806 else if (numBonds == 2)
2808 else if (numBonds == 3)
2814 strcmp(atomType,
"NY" )==0 ||
2815 strcmp(atomType,
"NC2" )==0 ||
2816 strcmp(atomType,
"N" )==0 ||
2817 strcmp(atomType,
"NH1" )==0 ||
2818 strcmp(atomType,
"NH2" )==0 ||
2819 strcmp(atomType,
"NR1" )==0 ||
2820 strcmp(atomType,
"NR2" )==0 ||
2821 strcmp(atomType,
"NR3" )==0 ||
2822 strcmp(atomType,
"NPH" )==0 ||
2823 strcmp(atomType,
"NC" )==0
2827 else if (numBonds == 2)
2829 else if (numBonds == 3)
2838 }
else if (atomType[0] ==
'O') {
2840 strcmp(atomType,
"OH1" )==0 ||
2841 strcmp(atomType,
"OS" )==0 ||
2842 strcmp(atomType,
"OC" )==0 ||
2843 strcmp(atomType,
"OT" )==0
2847 else if (numBonds == 2)
2852 strcmp(atomType,
"O" )==0 ||
2853 strcmp(atomType,
"OB" )==0 ||
2854 strcmp(atomType,
"OST" )==0 ||
2855 strcmp(atomType,
"OCA" )==0 ||
2856 strcmp(atomType,
"OM" )==0
2860 strcmp(atomType,
"OC" )==0
2868 }
else if (atomType[0] ==
'S') {
2875 }
else if (atomType[0] ==
'P') {
2878 else if (numBonds == 4)
2893 void Molecule::assignLCPOTypes(
int inputType) {
2894 int *heavyBonds =
new int[numAtoms];
2895 for (
int i = 0; i < numAtoms; i++)
2897 for (
int i = 0; i < numBonds; i++ ) {
2898 Bond curBond = bonds[i];
2899 int a1 = bonds[i].
atom1;
2900 int a2 = bonds[i].atom2;
2901 if (
atoms[a1].mass > 2.f &&
atoms[a2].mass > 2.f) {
2907 lcpoParamType =
new int[numAtoms];
2910 for (
int i = 0; i < numAtoms; i++) {
2913 if (inputType == 1) {
2927 if (
atoms[i].mass < 1.5 && lcpoParamType[i] != 0 ) {
2930 lcpoParamType[i] = 0;
2934 lcpoParamType[i] = 0;
2937 CkPrintf(
"ERROR in Molecule::assignLCPOTypes(): "
2938 "Light atom given heavy atom LCPO type.\n");
2946 iout <<
iWARN <<
"LONE PAIRS TO BE IGNORED BY SASA\n" <<
endi;
2949 iout <<
iWARN <<
"DRUDE PARTICLES TO BE IGNORED BY SASA\n" <<
endi;
2952 delete [] heavyBonds;
2956 void Molecule::plgLoadAtomBasics(molfile_atom_t *atomarray){
2962 hydrogenGroup.resize(0);
2966 for(
int i=0; i<numAtoms; i++) {
2967 int reslength = strlen(atomarray[i].resname)+1;
2968 int namelength = strlen(atomarray[i].name)+1;
2969 int typelength = strlen(atomarray[i].type)+1;
2970 atomNames[i].resname = nameArena->getNewArray(reslength);
2971 atomNames[i].atomname = nameArena->getNewArray(namelength);
2972 atomNames[i].atomtype = nameArena->getNewArray(typelength);
2973 strcpy(atomNames[i].resname, atomarray[i].resname);
2974 strcpy(atomNames[i].atomname, atomarray[i].name);
2975 strcpy(atomNames[i].atomtype, atomarray[i].type);
2977 atoms[i].mass = atomarray[i].mass;
2978 atoms[i].charge = atomarray[i].charge;
2983 tmpResLookup = tmpResLookup->
append(atomarray[i].segid, atomarray[i].resid, i);
2988 memcpy(one->
segname, atomarray[i].segid, strlen(atomarray[i].segid)+1);
2989 one->
resid = atomarray[i].resid;
2993 }
else if(
atoms[i].mass <= 0.05) {
2995 }
else if(
atoms[i].mass < 1.0) {
2997 }
else if(
atoms[i].mass <= 3.5) {
2999 }
else if((atomNames[i].atomname[0] ==
'O') &&
3000 (
atoms[i].mass>=14.0) && (
atoms[i].mass<=18.0)){
3008 void Molecule::plgLoadBonds(
int *from,
int *to){
3009 bonds =
new Bond[numBonds];
3010 int realNumBonds = 0;
3011 for(
int i=0; i<numBonds; i++) {
3012 Bond *thisBond = bonds+realNumBonds;
3013 thisBond->
atom1 = from[i]-1;
3014 thisBond->
atom2 = to[i]-1;
3017 atomNames[thisBond->
atom1].atomtype,
3018 atomNames[thisBond->
atom2].atomtype,
3024 if (is_lonepairs_psf) {
3026 if(k!=0. || is_lp(thisBond->
atom1) ||
3027 is_lp(thisBond->
atom2)) {
3031 if(k != 0.) realNumBonds++;
3035 if(numBonds != realNumBonds) {
3036 iout <<
iWARN <<
"Ignored" << numBonds-realNumBonds <<
3037 "bonds with zero force constants.\n" <<
endi;
3038 iout <<
iWARN <<
"Will get H-H distance in rigid H20 from H-O-H angle.\n" <<
endi;
3040 numBonds = realNumBonds;
3043 void Molecule::plgLoadAngles(
int *plgAngles)
3045 angles=
new Angle[numAngles];
3046 int *atomid = plgAngles;
3047 int numRealAngles = 0;
3048 for(
int i=0; i<numAngles; i++) {
3049 Angle *thisAngle = angles+numRealAngles;
3050 thisAngle->
atom1 = atomid[0]-1;
3051 thisAngle->
atom2 = atomid[1]-1;
3052 thisAngle->
atom3 = atomid[2]-1;
3056 atomNames[thisAngle->
atom1].atomtype,
3057 atomNames[thisAngle->
atom2].atomtype,
3058 atomNames[thisAngle->
atom3].atomtype,
3059 thisAngle, simParams->
alchOn ? -1 : 0);
3061 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n"
3065 Real k, t0, k_ub, r_ub;
3066 if ( thisAngle->
angle_type == -1 ) { k = -1.; k_ub = -1.; }
else
3068 if(k!=0. || k_ub!=0.) numRealAngles++;
3071 if(numAngles != numRealAngles) {
3072 iout <<
iWARN <<
"Ignored" << numAngles-numRealAngles <<
3073 " angles with zero force constants.\n" <<
endi;
3075 numAngles = numRealAngles;
3078 void Molecule::plgLoadDihedrals(
int *plgDihedrals)
3080 std::map< std::string, int > cache;
3083 int multiplicity = 1;
3085 lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
3086 dihedrals =
new Dihedral[numDihedrals];
3087 int numRealDihedrals = 0;
3088 int *atomid = plgDihedrals;
3089 for(
int i=0; i<numDihedrals; i++, atomid+=4) {
3090 Dihedral *thisDihedral = dihedrals + numRealDihedrals;
3092 for(
int j=0; j<4; j++) {
3093 if(atomid[j] != lastAtomIds[j]) {
3094 duplicate_bond =
FALSE;
3096 lastAtomIds[j] = atomid[j];
3099 if(duplicate_bond) {
3101 if(multiplicity==2) {
3102 numMultipleDihedrals++;
3109 thisDihedral->
atom1 = atomid[0]-1;
3110 thisDihedral->
atom2 = atomid[1]-1;
3111 thisDihedral->
atom3 = atomid[2]-1;
3112 thisDihedral->
atom4 = atomid[3]-1;
3115 sprintf(query,
"%10s :: %10s :: %10s :: %10s :: %d",
3116 atomNames[atomid[0]-1].atomtype,
3117 atomNames[atomid[1]-1].atomtype,
3118 atomNames[atomid[2]-1].atomtype,
3119 atomNames[atomid[3]-1].atomtype,
3121 auto search = cache.find(query);
3122 if ( search != cache.end() ) {
3126 atomNames[atomid[0]-1].atomtype,
3127 atomNames[atomid[1]-1].atomtype,
3128 atomNames[atomid[2]-1].atomtype,
3129 atomNames[atomid[3]-1].atomtype,
3130 thisDihedral, multiplicity, simParams->
alchOn ? -1 : 0);
3132 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n"
3139 numDihedrals = numRealDihedrals;
3142 void Molecule::plgLoadImpropers(
int *plgImpropers)
3145 int multiplicity = 1;
3147 lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
3148 impropers =
new Improper[numImpropers];
3149 int numRealImpropers = 0;
3150 int *atomid = plgImpropers;
3151 for(
int i=0; i<numImpropers; i++, atomid+=4) {
3152 Improper *thisImproper = impropers + numRealImpropers;
3154 for(
int j=0; j<4; j++) {
3155 if(atomid[j] != lastAtomIds[j]) {
3156 duplicate_bond =
FALSE;
3158 lastAtomIds[j] = atomid[j];
3161 if(duplicate_bond) {
3163 if(multiplicity==2) {
3164 numMultipleImpropers++;
3171 thisImproper->
atom1 = atomid[0]-1;
3172 thisImproper->
atom2 = atomid[1]-1;
3173 thisImproper->
atom3 = atomid[2]-1;
3174 thisImproper->
atom4 = atomid[3]-1;
3177 atomNames[atomid[0]-1].atomtype,
3178 atomNames[atomid[1]-1].atomtype,
3179 atomNames[atomid[2]-1].atomtype,
3180 atomNames[atomid[3]-1].atomtype,
3181 thisImproper, multiplicity);
3184 numImpropers = numRealImpropers;
3187 void Molecule::plgLoadCrossterms(
int *plgCterms)
3191 for(
int i=0; i<8; i++)
3194 crossterms =
new Crossterm[numCrossterms];
3195 int numRealCrossterms = 0;
3196 int *atomid = plgCterms;
3197 for(
int i=0; i<numCrossterms; i++, atomid+=8) {
3198 Crossterm *thisCrossterm = crossterms + numRealCrossterms;
3200 for(
int j=0; j<8; j++) {
3201 if(atomid[j] != lastAtomIds[j]) {
3202 duplicate_bond =
FALSE;
3204 lastAtomIds[j] = atomid[j];
3207 if(duplicate_bond) {
3210 numRealCrossterms++;
3212 thisCrossterm->
atom1 = atomid[0]-1;
3213 thisCrossterm->
atom2 = atomid[1]-1;
3214 thisCrossterm->
atom3 = atomid[2]-1;
3215 thisCrossterm->
atom4 = atomid[3]-1;
3216 thisCrossterm->
atom5 = atomid[4]-1;
3217 thisCrossterm->
atom6 = atomid[5]-1;
3218 thisCrossterm->
atom7 = atomid[6]-1;
3219 thisCrossterm->
atom8 = atomid[7]-1;
3222 atomNames[atomid[0]-1].atomtype,
3223 atomNames[atomid[1]-1].atomtype,
3224 atomNames[atomid[2]-1].atomtype,
3225 atomNames[atomid[3]-1].atomtype,
3226 atomNames[atomid[4]-1].atomtype,
3227 atomNames[atomid[5]-1].atomtype,
3228 atomNames[atomid[6]-1].atomtype,
3229 atomNames[atomid[7]-1].atomtype,
3233 numCrossterms = numRealCrossterms;
3237 occupancy =
new float[numAtoms];
3238 for(
int i=0; i<numAtoms; i++) {
3239 occupancy[i] = atomarray[i].occupancy;
3244 bfactor =
new float[numAtoms];
3245 for(
int i=0; i<numAtoms; i++) {
3246 bfactor[i] = atomarray[i].bfactor;
3261 void Molecule::build_lists_by_atom()
3265 register int numFixedAtoms = this->numFixedAtoms;
3275 bondsWithAtom =
new int32 *[numAtoms];
3276 cluster =
new int32 [numAtoms];
3277 clusterSize =
new int32 [numAtoms];
3279 bondsByAtom =
new int32 *[numAtoms];
3280 anglesByAtom =
new int32 *[numAtoms];
3281 dihedralsByAtom =
new int32 *[numAtoms];
3282 impropersByAtom =
new int32 *[numAtoms];
3283 crosstermsByAtom =
new int32 *[numAtoms];
3286 fullExclusionsByAtom =
new int32 *[numAtoms];
3287 modExclusionsByAtom =
new int32 *[numAtoms];
3290 gromacsPairByAtom =
new int32 *[numAtoms];
3295 const int pair_self =
3298 DebugM(3,
"Building bond lists.\n");
3301 for (i=0; i<numAtoms; i++)
3305 for (i=0; i<numRealBonds; i++)
3307 byAtomSize[bonds[i].atom1]++;
3308 byAtomSize[bonds[i].atom2]++;
3310 for (i=0; i<numAtoms; i++)
3312 bondsWithAtom[i] = tmpArena->getNewArray(byAtomSize[i]+1);
3313 bondsWithAtom[i][byAtomSize[i]] = -1;
3316 for (i=0; i<numRealBonds; i++)
3318 int a1 = bonds[i].atom1;
3319 int a2 = bonds[i].atom2;
3320 bondsWithAtom[a1][byAtomSize[a1]++] = i;
3321 bondsWithAtom[a2][byAtomSize[a2]++] = i;
3329 DebugM(3,
"Calculating exclusions for QM simulation.\n");
3332 delete_qm_bonded() ;
3334 DebugM(3,
"Re-Building bond lists.\n");
3338 for (i=0; i<numAtoms; i++)
3342 for (i=0; i<numRealBonds; i++)
3344 byAtomSize[bonds[i].atom1]++;
3345 byAtomSize[bonds[i].atom2]++;
3347 for (i=0; i<numAtoms; i++)
3349 bondsWithAtom[i][byAtomSize[i]] = -1;
3352 for (i=0; i<numRealBonds; i++)
3354 int a1 = bonds[i].atom1;
3355 int a2 = bonds[i].atom2;
3356 bondsWithAtom[a1][byAtomSize[a1]++] = i;
3357 bondsWithAtom[a2][byAtomSize[a2]++] = i;
3362 for (i=0; i<numAtoms; i++) {
3365 for (i=0; i<numAtoms; i++) {
3367 while ( cluster[ci] != ci ) ci = cluster[ci];
3368 for (
int32 *b = bondsWithAtom[i]; *b != -1; ++b ) {
3369 int a = bonds[*b].atom1;
3370 if ( a == i ) a = bonds[*b].atom2;
3373 while ( cluster[ca] != ca ) ca = cluster[ca];
3374 if ( ca > ci ) cluster[ca] = cluster[ci];
3375 else cluster[ci] = cluster[ca];
3381 for (i=0; i<numAtoms; i++) {
3382 int ci = cluster[i];
3383 if ( cluster[ci] != ci ) {
3385 cluster[i] = cluster[ci];
3391 for (i=0; i<numAtoms; i++) {
3394 for (i=0; i<numAtoms; i++) {
3395 clusterSize[cluster[i]] += 1;
3408 for (i=0; i<numAtoms; i++)
3413 for (i=0; i<numBonds; i++)
3415 if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3416 && fixedAtomFlags[bonds[i].atom2] )
continue;
3418 if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1)
continue;
3419 byAtomSize[bonds[i].atom1]++;
3422 for (i=0; i<numAtoms; i++)
3424 bondsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3425 bondsByAtom[i][byAtomSize[i]] = -1;
3428 for (i=0; i<numBonds; i++)
3430 if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3431 && fixedAtomFlags[bonds[i].atom2] )
continue;
3432 if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1)
continue;
3433 int a1 = bonds[i].atom1;
3434 bondsByAtom[a1][byAtomSize[a1]++] = i;
3436 for (i=0; i<numBonds; ++i) {
3437 int a1 = bonds[i].atom1;
3438 int a2 = bonds[i].atom2;
3442 sprintf(buff,
"Atom %d is bonded to itself", a1+1);
3445 for ( j = 0; j < byAtomSize[a1]; ++j ) {
3446 int b = bondsByAtom[a1][j];
3447 int ba1 = bonds[b].atom1;
3448 int ba2 = bonds[b].atom2;
3449 if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3451 sprintf(buff,
"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3455 for ( j = 0; j < byAtomSize[a2]; ++j ) {
3456 int b = bondsByAtom[a2][j];
3457 int ba1 = bonds[b].atom1;
3458 int ba2 = bonds[b].atom2;
3459 if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3461 sprintf(buff,
"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3467 DebugM(3,
"Building angle lists.\n");
3470 for (i=0; i<numAtoms; i++)
3475 for (i=0; i<numAngles; i++)
3477 if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3478 && fixedAtomFlags[angles[i].atom2]
3479 && fixedAtomFlags[angles[i].atom3] )
continue;
3480 if ( pair_self && fepAtomFlags[angles[i].atom1] != 1)
continue;
3481 byAtomSize[angles[i].atom1]++;
3484 for (i=0; i<numAtoms; i++)
3486 anglesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3487 anglesByAtom[i][byAtomSize[i]] = -1;
3490 for (i=0; i<numAngles; i++)
3492 if ( pair_self && fepAtomFlags[angles[i].atom1] != 1)
continue;
3493 if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3494 && fixedAtomFlags[angles[i].atom2]
3495 && fixedAtomFlags[angles[i].atom3] )
continue;
3496 int a1 = angles[i].atom1;
3497 anglesByAtom[a1][byAtomSize[a1]++] = i;
3500 DebugM(3,
"Building improper lists.\n");
3503 for (i=0; i<numAtoms; i++)
3507 numCalcImpropers = 0;
3508 for (i=0; i<numImpropers; i++)
3510 if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3511 && fixedAtomFlags[impropers[i].atom2]
3512 && fixedAtomFlags[impropers[i].atom3]
3513 && fixedAtomFlags[impropers[i].atom4] )
continue;
3514 if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1)
continue;
3515 byAtomSize[impropers[i].atom1]++;
3518 for (i=0; i<numAtoms; i++)
3520 impropersByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3521 impropersByAtom[i][byAtomSize[i]] = -1;
3524 for (i=0; i<numImpropers; i++)
3526 if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3527 && fixedAtomFlags[impropers[i].atom2]
3528 && fixedAtomFlags[impropers[i].atom3]
3529 && fixedAtomFlags[impropers[i].atom4] )
continue;
3530 if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1)
continue;
3531 int a1 = impropers[i].atom1;
3532 impropersByAtom[a1][byAtomSize[a1]++] = i;
3535 DebugM(3,
"Building dihedral lists.\n");
3538 for (i=0; i<numAtoms; i++)
3542 numCalcDihedrals = 0;
3543 for (i=0; i<numDihedrals; i++)
3545 if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3546 && fixedAtomFlags[dihedrals[i].atom2]
3547 && fixedAtomFlags[dihedrals[i].atom3]
3548 && fixedAtomFlags[dihedrals[i].atom4] )
continue;
3549 if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1)
continue;
3550 byAtomSize[dihedrals[i].atom1]++;
3553 for (i=0; i<numAtoms; i++)
3555 dihedralsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3556 dihedralsByAtom[i][byAtomSize[i]] = -1;
3559 for (i=0; i<numDihedrals; i++)
3561 if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3562 && fixedAtomFlags[dihedrals[i].atom2]
3563 && fixedAtomFlags[dihedrals[i].atom3]
3564 && fixedAtomFlags[dihedrals[i].atom4] )
continue;
3565 if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1)
continue;
3566 int a1 = dihedrals[i].atom1;
3567 dihedralsByAtom[a1][byAtomSize[a1]++] = i;
3570 DebugM(3,
"Building crossterm lists.\n");
3573 for (i=0; i<numAtoms; i++)
3577 numCalcCrossterms = 0;
3578 for (i=0; i<numCrossterms; i++)
3580 if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3581 && fixedAtomFlags[crossterms[i].atom2]
3582 && fixedAtomFlags[crossterms[i].atom3]
3583 && fixedAtomFlags[crossterms[i].atom4]
3584 && fixedAtomFlags[crossterms[i].atom5]
3585 && fixedAtomFlags[crossterms[i].atom6]
3586 && fixedAtomFlags[crossterms[i].atom7]
3587 && fixedAtomFlags[crossterms[i].atom8] )
continue;
3588 if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1)
continue;
3589 byAtomSize[crossterms[i].atom1]++;
3590 numCalcCrossterms++;
3592 for (i=0; i<numAtoms; i++)
3594 crosstermsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3595 crosstermsByAtom[i][byAtomSize[i]] = -1;
3598 for (i=0; i<numCrossterms; i++)
3600 if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3601 && fixedAtomFlags[crossterms[i].atom2]
3602 && fixedAtomFlags[crossterms[i].atom3]
3603 && fixedAtomFlags[crossterms[i].atom4]
3604 && fixedAtomFlags[crossterms[i].atom5]
3605 && fixedAtomFlags[crossterms[i].atom6]
3606 && fixedAtomFlags[crossterms[i].atom7]
3607 && fixedAtomFlags[crossterms[i].atom8] )
continue;
3608 if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1)
continue;
3609 int a1 = crossterms[i].atom1;
3610 crosstermsByAtom[a1][byAtomSize[a1]++] = i;
3614 if (is_lonepairs_psf) {
3616 DebugM(3,
"Initializing lone pair host index array.\n");
3617 lphostIndexes =
new int32[numAtoms];
3618 for (i = 0; i < numAtoms; i++) {
3619 lphostIndexes[i] = -1;
3621 for (i = 0; i < numLphosts; i++) {
3622 int32 index = lphosts[i].atom1;
3623 lphostIndexes[index] = i;
3629 DebugM(3,
"Building gromacsPair lists.\n");
3632 for (i=0; i<numAtoms; i++)
3639 if ( numFixedAtoms && fixedAtomFlags[
gromacsPair[i].atom1]
3640 && fixedAtomFlags[
gromacsPair[i].atom2] )
continue;
3641 if ( pair_self && fepAtomFlags[
gromacsPair[i].atom1] != 1)
continue;
3645 for (i=0; i<numAtoms; i++)
3647 gromacsPairByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3648 gromacsPairByAtom[i][byAtomSize[i]] = -1;
3653 if ( numFixedAtoms && fixedAtomFlags[
gromacsPair[i].atom1]
3654 && fixedAtomFlags[
gromacsPair[i].atom2] )
continue;
3655 if ( pair_self && fepAtomFlags[
gromacsPair[i].atom1] != 1)
continue;
3657 gromacsPairByAtom[a1][byAtomSize[a1]++] = i;
3662 DebugM(3,
"Building exclusion data.\n");
3669 delete [] bondsWithAtom; bondsWithAtom = 0;
3670 delete tmpArena; tmpArena = 0;
3676 numTotalExclusions = exclusionSet.size();
3678 iout <<
iINFO <<
"ADDED " << (numTotalExclusions - numExclusions) <<
" IMPLICIT EXCLUSIONS\n" <<
endi;
3682 for ( exclIter=exclIter.begin(),i=0; exclIter != exclIter.end(); exclIter++,i++ )
3688 exclusionSet.clear();
3690 DebugM(3,
"Building exclusion lists.\n");
3692 for (i=0; i<numAtoms; i++)
3696 numCalcExclusions = 0;
3697 numCalcFullExclusions = 0;
3698 for (i=0; i<numTotalExclusions; i++)
3700 if ( numFixedAtoms && fixedAtomFlags[
exclusions[i].atom1]
3701 && fixedAtomFlags[
exclusions[i].atom2] )
continue;
3703 numCalcExclusions++;
3704 if ( ! exclusions[i].modified ) numCalcFullExclusions++;
3707 for (i=0; i<numAtoms; i++)
3713 for (i=0; i<numTotalExclusions; i++)
3715 if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3716 && fixedAtomFlags[exclusions[i].atom2] )
continue;
3717 int a1 = exclusions[i].atom1;
3723 for (i=0; i<numAtoms; i++)
3729 for (i=0; i<numTotalExclusions; i++)
3731 if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3732 && fixedAtomFlags[exclusions[i].atom2] )
continue;
3733 if ( exclusions[i].modified ) {
3734 byAtomSize2[exclusions[i].atom1]++;
3735 byAtomSize2[exclusions[i].atom2]++;
3737 byAtomSize[exclusions[i].atom1]++;
3738 byAtomSize[exclusions[i].atom2]++;
3742 for (i=0; i<numAtoms; i++)
3744 fullExclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3745 fullExclusionsByAtom[i][0] = 0;
3746 modExclusionsByAtom[i] = arena->getNewArray(byAtomSize2[i]+1);
3747 modExclusionsByAtom[i][0] = 0;
3750 for (i=0; i<numTotalExclusions; i++)
3752 int a1 = exclusions[i].atom1;
3753 int a2 = exclusions[i].atom2;
3754 if ( numFixedAtoms && fixedAtomFlags[a1]
3755 && fixedAtomFlags[a2] )
continue;
3757 if ( exclusions[i].modified ) {
3758 l1 = modExclusionsByAtom[a1];
3759 l2 = modExclusionsByAtom[a2];
3761 l1 = fullExclusionsByAtom[a1];
3762 l2 = fullExclusionsByAtom[a2];
3769 for (i=0; i<numAtoms; i++) {
3770 int32 *lf = fullExclusionsByAtom[i];
3771 iout <<
"EXCL " << i <<
" FULL";
3773 for (
int j = 0; j < nf; ++j ) {
3774 iout <<
" " << *(lf++);
3777 int32 *lm = modExclusionsByAtom[i];
3778 iout <<
"EXCL " << i <<
" MOD";
3780 for (
int j = 0; j < nm; ++j ) {
3781 iout <<
" " << *(lm++);
3788 if (is_drude_psf || simParams->
drudeOn) {
3794 if (tholes != NULL)
delete[] tholes;
3798 for (i = 0; i < numTotalExclusions; i++) {
3800 if (exclusions[i].modified)
continue;
3801 int a1 = exclusions[i].atom1;
3802 int a2 = exclusions[i].atom2;
3803 if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3809 if (numTholes != 0) tholes =
new Thole[numTholes];
3816 for (i = 0; i < numTotalExclusions; i++) {
3818 if (exclusions[i].modified)
continue;
3819 int a1 = exclusions[i].atom1;
3820 int a2 = exclusions[i].atom2;
3822 if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3823 Real thsum = drudeConsts[a1].thole + drudeConsts[a2].thole;
3824 Real aprod = drudeConsts[a1].alpha * drudeConsts[a2].alpha;
3826 Real apower = (aprod <= 0 ? 0 : powf(aprod, -1.f/6));
3827 tholes[nt].atom1 = a1;
3828 tholes[nt].atom2 = a1+1;
3829 tholes[nt].atom3 = a2;
3830 tholes[nt].atom4 = a2+1;
3831 tholes[nt].aa = apower * thsum;
3832 tholes[nt].qq = c *
atoms[a1+1].charge *
atoms[a2+1].charge;
3838 DebugM(3,
"Building Thole correction term lists.\n");
3839 tholesByAtom =
new int32 *[numAtoms];
3841 for (i = 0; i < numAtoms; i++) {
3845 for (i = 0; i < numTholes; i++) {
3846 if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
3847 && fixedAtomFlags[tholes[i].atom2]
3848 && fixedAtomFlags[tholes[i].atom3]
3849 && fixedAtomFlags[tholes[i].atom4] )
continue;
3850 if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1)
continue;
3851 byAtomSize[tholes[i].atom1]++;
3854 for (i = 0; i < numAtoms; i++) {
3855 tholesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3856 tholesByAtom[i][byAtomSize[i]] = -1;
3859 for (i = 0; i < numTholes; i++) {
3860 if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
3861 && fixedAtomFlags[tholes[i].atom2]
3862 && fixedAtomFlags[tholes[i].atom3]
3863 && fixedAtomFlags[tholes[i].atom4] )
continue;
3864 if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1)
continue;
3865 int a1 = tholes[i].atom1;
3866 tholesByAtom[a1][byAtomSize[a1]++] = i;
3870 DebugM(3,
"Building anisotropic term lists.\n");
3871 anisosByAtom =
new int32 *[numAtoms];
3873 for (i = 0; i < numAtoms; i++) {
3877 for (i = 0; i < numAnisos; i++) {
3878 if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
3879 && fixedAtomFlags[anisos[i].atom2]
3880 && fixedAtomFlags[anisos[i].atom3]
3881 && fixedAtomFlags[anisos[i].atom4] )
continue;
3882 if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1)
continue;
3883 byAtomSize[anisos[i].atom1]++;
3886 for (i = 0; i < numAtoms; i++) {
3887 anisosByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3888 anisosByAtom[i][byAtomSize[i]] = -1;
3891 for (i = 0; i < numAnisos; i++) {
3892 if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
3893 && fixedAtomFlags[anisos[i].atom2]
3894 && fixedAtomFlags[anisos[i].atom3]
3895 && fixedAtomFlags[anisos[i].atom4] )
continue;
3896 if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1)
continue;
3897 int a1 = anisos[i].atom1;
3898 anisosByAtom[a1][byAtomSize[a1]++] = i;
3904 delete [] byAtomSize; byAtomSize = 0;
3905 delete [] byAtomSize2; byAtomSize2 = 0;
3911 for (i=0; i<numAtoms; i++)
3913 all_exclusions[i].
min = numAtoms;
3914 all_exclusions[i].max = -1;
3916 for (i=0; i<numTotalExclusions; i++)
3919 int a1 = exclusions[i].atom1;
3920 int a2 = exclusions[i].atom2;
3921 if ( numFixedAtoms && fixedAtomFlags[a1]
3922 && fixedAtomFlags[a2] )
continue;
3923 if ( all_exclusions[a1].min > a2 ) all_exclusions[a1].min = a2;
3924 if ( all_exclusions[a2].min > a1 ) all_exclusions[a2].min = a1;
3925 if ( a2 > all_exclusions[a1].max ) all_exclusions[a1].max = a2;
3926 if ( a1 > all_exclusions[a2].max ) all_exclusions[a2].max = a1;
3930 int maxDenseAllFull = 0;
3931 int numDenseAllFull = 0;
3932 for (i=0; i<numAtoms; i++) {
3933 int iInMiddle = ( i < all_exclusions[i].max &&
3934 i > all_exclusions[i].min ) ? 1 : 0;
3935 int s = all_exclusions[i].max - all_exclusions[i].min + 1;
3936 if ( s == fullExclusionsByAtom[i][0] + iInMiddle ) {
3937 if ( s > maxDenseAllFull ) maxDenseAllFull = s;
3938 all_exclusions[i].flags = (
char*)-1;
3940 all_exclusions[i].flags = 0;
3943 char *denseFullArray = exclArena->getNewArray(maxDenseAllFull);
3944 for ( i=0; i<maxDenseAllFull; ++i ) denseFullArray[i] =
EXCHCK_FULL;
3946 int exclmem = maxDenseAllFull;
3948 for (i=0; i<numAtoms; i++) {
3949 int s = all_exclusions[i].max - all_exclusions[i].min + 1;
3950 if ( all_exclusions[i].max != -1 ) {
3951 if ( all_exclusions[i].flags ) {
3952 all_exclusions[i].flags = denseFullArray;
3954 }
else if ( s < maxExclusionFlags ) {
3955 char *f = all_exclusions[i].flags = exclArena->getNewArray(s);
3956 for (
int k=0; k<s; ++k ) f[k] = 0;
3959 all_exclusions[i].flags = 0;
3962 all_exclusions[i].flags = (
char*)-1;
3966 iout <<
iINFO << numTotalExclusions <<
" exclusions consume "
3967 << exclmem <<
" bytes.\n" <<
endi;
3969 <<
" atoms sharing one array.\n" <<
endi;
3971 for (i=0; i<numTotalExclusions; i++)
3973 int a1 = exclusions[i].atom1;
3974 int a2 = exclusions[i].atom2;
3975 if ( numFixedAtoms && fixedAtomFlags[a1]
3976 && fixedAtomFlags[a2] )
continue;
3977 if ( exclusions[i].modified ) {
3978 if ( all_exclusions[a1].flags )
3979 all_exclusions[a1].flags[a2-all_exclusions[a1].min] =
EXCHCK_MOD;
3980 if ( all_exclusions[a2].flags )
3981 all_exclusions[a2].flags[a1-all_exclusions[a2].min] =
EXCHCK_MOD;
3983 if ( all_exclusions[a1].flags )
3984 all_exclusions[a1].flags[a2-all_exclusions[a1].min] =
EXCHCK_FULL;
3985 if ( all_exclusions[a2].flags )
3986 all_exclusions[a2].flags[a1-all_exclusions[a2].min] =
EXCHCK_FULL;
4015 void Molecule::build_exclusions()
4020 exclude_flag = simParams->
exclude;
4023 for (i=0; i<numExclusions; i++)
4025 exclusionSet.add(exclusions[i]);
4033 switch (exclude_flag)
4060 if (is_lonepairs_psf || is_drude_psf) {
4061 build_inherited_excl(
SCALED14 == exclude_flag);
4080 void Molecule::build_inherited_excl(
int modified) {
4082 int32 *bond1, *bond2, *bond3, *bond4, *bond5;
4083 int32 i, j, mid1, mid2, mid3, mid4;
4087 for (i = 0; i < numAtoms; i++) {
4089 if (!is_drude(i) && !is_lp(i))
continue;
4093 bond1 = bondsWithAtom[i];
4097 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4098 sprintf(err_msg,
"FOUND ISOLATED %s PARTICLE %d", idescrip, i+1);
4101 if (-1 != *(bond1+1)) {
4103 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4104 sprintf(err_msg,
"FOUND MULTIPLY LINKED %s PARTICLE %d",
4110 mid1 = bonds[*bond1].atom1;
4111 if (mid1 == i) mid1 = bonds[*bond1].atom2;
4114 if (is_drude(mid1) || is_lp(mid1) || is_hydrogen(mid1)) {
4116 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4117 sprintf(err_msg,
"PARENT ATOM %d of %s PARTICLE %d "
4118 "IS NOT HEAVY ATOM", mid1+1, idescrip, i+1);
4122 if (exclude_flag ==
NONE) {
4132 bond2 = bondsWithAtom[mid1];
4133 while (*bond2 != -1) {
4134 j = bonds[*bond2].atom1;
4135 if ((is_drude(j) || is_lp(j)) && j != mid1) {
4136 if (i < j) exclusionSet.add(
Exclusion(i, j));
4137 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4139 j = bonds[*bond2].atom2;
4140 if ((is_drude(j) || is_lp(j)) && j != mid1) {
4141 if (i < j) exclusionSet.add(
Exclusion(i, j));
4142 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4150 bond2 = bondsWithAtom[mid1];
4153 while (*bond2 != -1) {
4154 if (bonds[*bond2].atom1 == mid1) {
4155 mid2 = bonds[*bond2].atom2;
4158 mid2 = bonds[*bond2].atom1;
4168 if (exclude_flag ==
ONETWO) {
4178 bond3 = bondsWithAtom[mid2];
4179 while (*bond3 != -1) {
4180 j = bonds[*bond3].atom1;
4181 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4182 if (i < j) exclusionSet.add(
Exclusion(i, j));
4183 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4185 j = bonds[*bond3].atom2;
4186 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4187 if (i < j) exclusionSet.add(
Exclusion(i, j));
4188 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4196 bond3 = bondsWithAtom[mid2];
4199 while (*bond3 != -1) {
4201 if (bonds[*bond3].atom1 == mid2) {
4202 mid3 = bonds[*bond3].atom2;
4205 mid3 = bonds[*bond3].atom1;
4219 else if (mid3 < i) {
4225 bond4 = bondsWithAtom[mid3];
4226 while (*bond4 != -1) {
4227 j = bonds[*bond4].atom1;
4228 if ((is_drude(j) || is_lp(j)) && j != mid3) {
4229 if (i < j) exclusionSet.add(
Exclusion(i, j));
4230 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4232 j = bonds[*bond4].atom2;
4233 if ((is_drude(j) || is_lp(j)) && j != mid3) {
4234 if (i < j) exclusionSet.add(
Exclusion(i, j));
4235 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4243 bond4 = bondsWithAtom[mid3];
4246 while (*bond4 != -1) {
4248 if (bonds[*bond4].atom1 == mid3) {
4249 mid4 = bonds[*bond4].atom2;
4252 mid4 = bonds[*bond4].atom1;
4262 if (is_drude(mid4) || is_lp(mid4)) {
4267 else if (mid4 < i) {
4278 int modi = modified;
4280 int amin = (mid1 < mid4 ? mid1 : mid4);
4281 int amax = (mid1 >= mid4 ? mid1 : mid4);
4293 exclusionSet.add(
Exclusion(i, mid4, modi));
4295 else if (mid4 < i) {
4296 exclusionSet.add(
Exclusion(mid4, i, modi));
4301 bond5 = bondsWithAtom[mid4];
4302 while (*bond5 != -1) {
4303 j = bonds[*bond5].atom1;
4304 if ((is_drude(j) || is_lp(j)) && j != mid4) {
4305 if (i<j) exclusionSet.add(
Exclusion(i,j,modi));
4306 else if (j<i) exclusionSet.add(
Exclusion(j,i,modi));
4308 j = bonds[*bond5].atom2;
4309 if ((is_drude(j) || is_lp(j)) && j != mid4) {
4310 if (i<j) exclusionSet.add(
Exclusion(i,j,modi));
4311 else if (j<i) exclusionSet.add(
Exclusion(j,i,modi));
4341 void Molecule::build12excl(
void)
4347 for (i=0; i<numAtoms; i++)
4349 current_val = bondsWithAtom[i];
4352 while (*current_val != -1)
4354 if (bonds[*current_val].atom1 == i)
4356 if (i<bonds[*current_val].atom2)
4358 exclusionSet.add(
Exclusion(i,bonds[*current_val].atom2));
4363 if (i<bonds[*current_val].atom1)
4365 exclusionSet.add(
Exclusion(i,bonds[*current_val].atom1));
4381 void Molecule::build13excl(
void)
4383 int32 *bond1, *bond2;
4389 for (i=0; i<numAtoms; i++)
4391 bond1 = bondsWithAtom[i];
4394 while (*bond1 != -1)
4396 if (bonds[*bond1].atom1 == i)
4398 middle_atom=bonds[*bond1].atom2;
4402 middle_atom=bonds[*bond1].atom1;
4405 bond2 = bondsWithAtom[middle_atom];
4409 while (*bond2 != -1)
4411 if (bonds[*bond2].atom1 == middle_atom)
4413 if (i < bonds[*bond2].atom2)
4415 exclusionSet.add(
Exclusion(i,bonds[*bond2].atom2));
4420 if (i < bonds[*bond2].atom1)
4422 exclusionSet.add(
Exclusion(i,bonds[*bond2].atom1));
4442 void Molecule::build14excl(
int modified)
4444 int32 *bond1, *bond2, *bond3;
4449 for (i=0; i<numAtoms; i++)
4451 if (is_drude(i) || is_lp(i))
continue;
4454 bond1 = bondsWithAtom[i];
4456 while (*bond1 != -1)
4458 if (bonds[*bond1].atom1 == i)
4460 mid1=bonds[*bond1].atom2;
4464 mid1=bonds[*bond1].atom1;
4467 bond2 = bondsWithAtom[mid1];
4470 while (*bond2 != -1)
4472 if (bonds[*bond2].atom1 == mid1)
4474 mid2 = bonds[*bond2].atom2;
4478 mid2 = bonds[*bond2].atom1;
4490 bond3=bondsWithAtom[mid2];
4493 while (*bond3 != -1)
4495 if (bonds[*bond3].atom1 == mid2)
4497 int j = bonds[*bond3].atom2;
4503 if (i < j && !is_drude(j) && !is_lp(j))
4505 exclusionSet.add(
Exclusion(i,j,modified));
4510 int j = bonds[*bond3].atom1;
4516 if (i < j && !is_drude(j) && !is_lp(j))
4518 exclusionSet.add(
Exclusion(i,j,modified));
4540 void Molecule::stripFepExcl(
void)
4546 for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
4548 int t1 = get_fep_type(exclIter->atom1);
4549 int t2 = get_fep_type(exclIter->atom2);
4550 if (t1 && t2 && t1 !=t2 && abs(t1-t2) != 2) {
4551 fepExclusionSet.
add(*exclIter);
4555 for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
4557 int ifep_type = get_fep_type(exclIter->atom1);
4558 int jfep_type = get_fep_type(exclIter->atom2);
4561 if (ifep_type != 1 || jfep_type != 1) {
4562 fepExclusionSet.
add(*exclIter);
4567 if (!(ifep_type == 1 && jfep_type == 2) &&
4568 !(ifep_type == 2 && jfep_type == 1)) {
4569 fepExclusionSet.
add(*exclIter);
4576 for ( fepIter=fepIter.begin(); fepIter != fepIter.end(); fepIter++ )
4578 exclusionSet.del(*fepIter);
4592 if ( (psf_file =
Fopen(fname,
"r")) == NULL)
4595 sprintf(err_msg,
"UNABLE TO OPEN THE COMPRESSED .psf FILE %s", fname);
4604 NAMD_die(
"The compressed psf file format is incorrect, please re-generate!\n");
4606 float psfVer = 0.0f;
4607 sscanf(buffer,
"FORMAT VERSION: %f\n", &psfVer);
4609 NAMD_die(
"The compressed psf file format is incorrect, please re-generate!\n");
4614 NAMD_die(
"UNABLE TO FIND NSEGMENTNAMES");
4615 sscanf(buffer,
"%d", &segNamePoolSize);
4617 if(segNamePoolSize!=0)
4619 for(
int i=0; i<segNamePoolSize; i++){
4621 sscanf(buffer,
"%s", strBuf);
4622 segNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4626 for(
int i=0; i<segNamePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4631 NAMD_die(
"UNABLE TO FIND NRESIDUENAMES");
4632 sscanf(buffer,
"%d", &resNamePoolSize);
4634 if(resNamePoolSize!=0)
4636 for(
int i=0; i<resNamePoolSize; i++){
4638 sscanf(buffer,
"%s", strBuf);
4639 resNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4643 for(
int i=0; i<resNamePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4648 NAMD_die(
"UNABLE TO FIND NATOMNAMES");
4649 sscanf(buffer,
"%d", &atomNamePoolSize);
4650 if(atomNamePoolSize!=0)
4652 for(
int i=0; i<atomNamePoolSize; i++){
4654 sscanf(buffer,
"%s", strBuf);
4655 atomNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4661 NAMD_die(
"UNABLE TO FIND NATOMTYPES");
4662 sscanf(buffer,
"%d", &atomTypePoolSize);
4664 if(atomTypePoolSize!=0)
4666 for(
int i=0; i<atomTypePoolSize; i++){
4668 sscanf(buffer,
"%s", strBuf);
4669 atomTypePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4673 for(
int i=0; i<atomTypePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4678 NAMD_die(
"UNABLE TO FIND NCHARGES");
4679 sscanf(buffer,
"%d", &chargePoolSize);
4680 if(chargePoolSize!=0)
4681 atomChargePool =
new Real[chargePoolSize];
4682 for(
int i=0; i<chargePoolSize; i++){
4684 sscanf(buffer,
"%f", atomChargePool+i);
4689 NAMD_die(
"UNABLE TO FIND NMASSES");
4690 sscanf(buffer,
"%d", &massPoolSize);
4692 atomMassPool =
new Real[massPoolSize];
4693 for(
int i=0; i<massPoolSize; i++){
4695 sscanf(buffer,
"%f", atomMassPool+i);
4700 NAMD_die(
"UNABLE TO FIND ATOMSIGS");
4701 sscanf(buffer,
"%d", &atomSigPoolSize);
4704 int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
4707 for(
int i=0; i<atomSigPoolSize; i++){
4711 NAMD_die(
"UNABLE TO FIND NBONDSIGS");
4712 sscanf(buffer,
"%d", &typeCnt);
4717 for(
int j=0; j<typeCnt; j++){
4719 sscanf(buffer,
"%d | %d | %d", &tmp1, &ttype, &tisReal);
4721 oneSig.offset[0] = tmp1;
4723 if(tisReal) numRealBonds++;
4729 NAMD_die(
"UNABLE TO FIND NTHETASIGS");
4730 sscanf(buffer,
"%d", &typeCnt);
4735 for(
int j=0; j<typeCnt; j++){
4737 sscanf(buffer,
"%d %d | %d | %d", &tmp1, &tmp2, &ttype, &tisReal);
4739 oneSig.offset[0] = tmp1;
4740 oneSig.offset[1] = tmp2;
4746 NAMD_die(
"UNABLE TO FIND NPHISIGS");
4747 sscanf(buffer,
"%d", &typeCnt);
4752 for(
int j=0; j<typeCnt; j++){
4754 sscanf(buffer,
"%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
4756 oneSig.offset[0] = tmp1;
4757 oneSig.offset[1] = tmp2;
4758 oneSig.offset[2] = tmp3;
4764 NAMD_die(
"UNABLE TO FIND NIMPHISIGS");
4765 sscanf(buffer,
"%d", &typeCnt);
4770 for(
int j=0; j<typeCnt; j++){
4772 sscanf(buffer,
"%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
4774 oneSig.offset[0] = tmp1;
4775 oneSig.offset[1] = tmp2;
4776 oneSig.offset[2] = tmp3;
4782 NAMD_die(
"UNABLE TO FIND NCRTERMSIGS");
4783 sscanf(buffer,
"%d", &typeCnt);
4788 for(
int j=0; j<typeCnt; j++){
4790 sscanf(buffer,
"%d %d %d %d %d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &ttype, &tisReal);
4792 oneSig.offset[0] = tmp1;
4793 oneSig.offset[1] = tmp2;
4794 oneSig.offset[2] = tmp3;
4795 oneSig.offset[3] = tmp4;
4796 oneSig.offset[4] = tmp5;
4797 oneSig.offset[5] = tmp6;
4798 oneSig.offset[6] = tmp7;
4806 NAMD_die(
"UNABLE TO FIND NEXCLSIGS");
4808 sscanf(buffer,
"%d", &exclSigPoolSize);
4810 vector<int> fullExcls;
4811 vector<int> modExcls;
4812 for(
int i=0; i<exclSigPoolSize; i++){
4814 for(
int j=0; j<fullExclCnt; j++)
4817 for(
int j=0; j<modExclCnt; j++)
4821 exclSigPool[i].setOffsets(fullExcls, modExcls);
4830 NAMD_die(
"UNABLE TO FIND NCLUSTERS");
4832 sscanf(buffer,
"%d", &numClusters);
4837 sscanf(buffer,
"%d", &numAtoms);
4841 NAMD_die(
"UNABLE TO FIND NHYDROGENGROUP");
4842 sscanf(buffer,
"%d", &numHydrogenGroups);
4846 NAMD_die(
"UNABLE TO FIND MAXHYDROGENGROUPSIZE");
4847 sscanf(buffer,
"%d", &maxHydrogenGroupSize);
4850 NAMD_die(
"UNABLE TO FIND NMIGRATIONGROUP");
4851 sscanf(buffer,
"%d", &numMigrationGroups);
4854 NAMD_die(
"UNABLE TO FIND MAXMIGRATIONGROUPSIZE");
4855 sscanf(buffer,
"%d", &maxMigrationGroupSize);
4857 int inputRigidType = -1;
4860 NAMD_die(
"UNABLE TO FIND RIGIDBONDTYPE");
4861 sscanf(buffer,
"%d", &inputRigidType);
4865 const char *tmpstr[]={
"RIGID_NONE",
"RIGID_ALL",
"RIGID_WATER"};
4867 sprintf(errmsg,
"RIGIDBOND TYPE MISMATCH BETWEEN INPUT (%s) AND CURRENT RUN (%s)",
4868 tmpstr[inputRigidType], tmpstr[simParams->
rigidBonds]);
4876 NAMD_die(
"UNABLE TO FIND OCCUPANCYVALID");
4877 sscanf(buffer,
"%d", &isOccupancyValid);
4880 NAMD_die(
"UNABLE TO FIND TEMPFACTORVALID");
4881 sscanf(buffer,
"%d", &isBFactorValid);
4889 build_extra_bonds(params, cfgList->
find(
"extraBondsFile"));
4893 NAMD_die(
"UNABLE TO FIND DIHEDRALPARAMARRAY");
4902 NAMD_die(
"UNABLE TO FIND IMPROPERPARAMARRAY");
4919 void Molecule::read_binary_atom_info(
int fromAtomID,
int toAtomID,
InputAtomList& inAtoms){
4920 int numAtomsPar = toAtomID-fromAtomID+1;
4921 CmiAssert(numAtomsPar > 0);
4922 CmiAssert(inAtoms.
size() == numAtomsPar);
4941 eachAtomMass = NULL;
4942 eachAtomCharge = NULL;
4944 eachAtomExclSig = NULL;
4967 FILE *perAtomFile = fopen(simParams->
binAtomFile,
"rb");
4968 if (perAtomFile==NULL) {
4970 sprintf(err_msg,
"UNABLE TO OPEN THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s", simParams->
binAtomFile);
4976 flipNum((
char *)&rMagicNum,
sizeof(
int), 1);
4978 fread(&fMagicNum,
sizeof(
int), 1, perAtomFile);
4979 if (fMagicNum==magicNum) {
4981 }
else if (fMagicNum==rMagicNum) {
4985 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS CORRUPTED", simParams->
binAtomFile);
4989 float verNum = 0.0f;
4990 fread(&verNum,
sizeof(
float), 1, perAtomFile);
4991 if (needFlip)
flipNum((
char *)&verNum,
sizeof(
float), 1);
4994 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n", simParams->
binAtomFile);
4999 fread(&recSize,
sizeof(
int), 1, perAtomFile);
5000 if(needFlip)
flipNum((
char *)&recSize,
sizeof(
int), 1);
5003 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM RECORD SIZE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n", simParams->
binAtomFile);
5007 const int BUFELEMS = 32*1024;
5012 if ( _fseeki64(perAtomFile,startbyte,SEEK_CUR) )
5014 if ( fseeko(perAtomFile,startbyte,SEEK_CUR) )
5018 sprintf(errmsg,
"Error on seeking binary file %s", simParams->
binAtomFile);
5024 int atomsCnt = numAtomsPar;
5027 while(atomsCnt >= BUFELEMS) {
5028 if ( fread((
char *)elemsBuf,
sizeof(
OutputAtomRecord), BUFELEMS, perAtomFile) != BUFELEMS ) {
5030 sprintf(errmsg,
"Error on reading binary file %s", simParams->
binAtomFile);
5034 for(
int i=0; i<BUFELEMS; i++, curIdx++, oneRec++) {
5036 int aid = curIdx+fromAtomID;
5037 if(needFlip) oneRec->
flip();
5038 load_one_inputatom(aid, oneRec, fAtom);
5040 atomsCnt -= BUFELEMS;
5043 if ( fread(elemsBuf,
sizeof(
OutputAtomRecord), atomsCnt, perAtomFile) != atomsCnt ) {
5045 sprintf(errmsg,
"Error on reading binary file %s", simParams->
binAtomFile);
5049 for(
int i=curIdx; i<numAtomsPar; i++, oneRec++) {
5051 int aid = i+fromAtomID;
5052 if(needFlip) oneRec->
flip();
5053 load_one_inputatom(aid,oneRec,fAtom);
5056 if ( fclose(perAtomFile) ) {
5058 sprintf(errmsg,
"Error on closing binary file %s", simParams->
binAtomFile);
5067 is_atom_fixed(fromAtomID, &listIdx);
5068 for(
int i=listIdx; i<fixedAtomsSet->size(); i++){
5069 const AtomSet one = fixedAtomsSet->item(i);
5071 int sAtomId = one.aid1>fromAtomID ? one.aid1:fromAtomID;
5072 int eAtomId = one.aid2>toAtomID? toAtomID:one.aid2;
5073 for(
int j=sAtomId; j<=eAtomId; j++)
5074 inAtoms[j-fromAtomID].atomFixed = 1;
5081 char *thisAtomName = NULL;
5143 }
else if (thisAtomMass <= 0.05) {
5145 }
else if (thisAtomMass < 1.0) {
5147 }
else if (thisAtomMass <= 3.5) {
5149 }
else if (thisAtomName[0]==
'O' &&
5150 (thisAtomMass >= 14.0) && (thisAtomMass <= 18.0)) {
5164 void Molecule::build_excl_check_signatures(){
5166 for(
int i=0; i<exclSigPoolSize; i++){
5174 int fullMin, fullMax, modMin, modMax;
5182 if(fullMin < modMin)
5183 sigChk->
min = fullMin;
5185 sigChk->
min = modMin;
5186 if(fullMax < modMax)
5187 sigChk->
max = modMax;
5189 sigChk->
max = fullMax;
5197 iout <<
iWARN <<
"an empty exclusion signature with index "
5198 << i <<
"!\n" <<
endi;
5203 sigChk->
flags =
new char[sigChk->
max-sigChk->
min+1];
5204 memset(sigChk->
flags, 0,
sizeof(
char)*(sigChk->
max-sigChk->
min+1));
5225 void Molecule::load_atom_set(
StringList *setfile,
const char *setname,
5226 int *numAtomsInSet, AtomSetList **atomsSet)
const {
5227 if(setfile == NULL) {
5229 sprintf(errmsg,
"The text input file for %s atoms is not found!", setname);
5232 FILE *ifp = fopen(setfile->
data,
"r");
5236 sprintf(errmsg,
"ERROR IN OPENING %s ATOMS FILE: %s\n", setname, setfile->
data);
5241 int numLocalAtoms = 0;
5243 AtomSetList *localAtomsSet =
new AtomSetList();
5248 bool hasDash =
false;
5249 for(
int i=0; oneline[i] && i<128; i++){
5250 if(oneline[i]==
'-') {
5256 sscanf(oneline,
"%d-%d", &(one.aid1), &(one.aid2));
5257 if(one.aid1>one.aid2 || one.aid1<0 || one.aid2<0) {
5259 sprintf(errmsg,
"The input for %s atoms is wrong: %s\n", setname, oneline);
5262 numLocalAtoms += (one.aid2-one.aid1+1);
5264 sscanf(oneline,
"%d", &(one.aid1));
5267 sprintf(errmsg,
"The input for %s atoms is wrong: %s\n", setname, oneline);
5270 one.aid2 = one.aid1;
5273 localAtomsSet->add(one);
5277 std::sort(localAtomsSet->begin(), localAtomsSet->end());
5279 *numAtomsInSet = numLocalAtoms;
5280 *atomsSet = localAtomsSet;
5283 void Molecule::load_fixed_atoms(
StringList *fixedfile){
5284 load_atom_set(fixedfile,
"FIXED", &numFixedAtoms, &fixedAtomsSet);
5287 void Molecule::load_constrained_atoms(
StringList *constrainedfile){
5288 load_atom_set(constrainedfile,
"CONSTRAINED", &numConstraints, &constrainedAtomsSet);
5291 Bool Molecule::is_atom_in_set(AtomSetList *localAtomsSet,
int aid,
int *listIdx)
const {
5292 int idx = localAtomsSet->size();
5294 int lIdx = localAtomsSet->size()-1;
5296 while(rIdx <= lIdx){
5297 int mIdx = (rIdx+lIdx)/2;
5298 const AtomSet one = localAtomsSet->item(mIdx);
5304 }
else if(aid > one.aid1){
5308 if(listIdx) *listIdx = mIdx;
5315 if(listIdx) *listIdx = mIdx;
5321 if(listIdx) *listIdx = idx;
5339 #ifdef MEM_OPT_VERSION
5340 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5349 <<
"******************************************\n" \
5350 <<
"NUM NAME TYPE RES MASS CHARGE CHARGE FEP-CHARGE" \
5351 <<
"SIGMA EPSILON SIGMA14 EPSILON14\n" \
5354 for (i=0; i<numAtoms; i++)
5359 DebugM(2,i+1 <<
" " << atomNames[i].atomname \
5360 <<
" " << atomNames[i].atomtype <<
" " \
5361 << atomNames[i].resname <<
" " << atoms[i].mass \
5362 <<
" " << atoms[i].charge <<
" " << sigma \
5363 <<
" " << epsilon <<
" " << sigma14 \
5364 <<
" " << epsilon14 <<
"\n" \
5382 #ifdef MEM_OPT_VERSION
5383 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5389 DebugM(2,
"BOND LIST\n" <<
"********************************\n" \
5390 <<
"ATOM1 ATOM2 TYPE1 TYPE2 k x0" \
5393 for (i=0; i<numBonds; i++)
5397 DebugM(2,bonds[i].atom1+1 <<
" " \
5398 << bonds[i].atom2+1 <<
" " \
5399 << atomNames[bonds[i].atom1].atomtype <<
" " \
5400 << atomNames[bonds[i].atom2].atomtype <<
" " << k \
5401 <<
" " << x0 <<
endi);
5419 #ifdef MEM_OPT_VERSION
5420 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5424 DebugM(2,
"EXPLICIT EXCLUSION LIST\n" \
5425 <<
"********************************\n" \
5429 for (i=0; i<numExclusions; i++)
5431 DebugM(2,exclusions[i].atom1+1 <<
" " \
5432 << exclusions[i].atom2+1 <<
endi);
5449 #ifdef MEM_OPT_VERSION
5455 msg->
put(massPoolSize);
5456 msg->
put(massPoolSize, atomMassPool);
5458 msg->
put(chargePoolSize);
5459 msg->
put(chargePoolSize, atomChargePool);
5462 msg->
put(atomSigPoolSize);
5463 for(
int i=0; i<atomSigPoolSize; i++)
5467 msg->
put(exclSigPoolSize);
5468 for(
int i=0; i<exclSigPoolSize; i++)
5469 exclSigPool[i].pack(msg);
5471 msg->
put(numHydrogenGroups);
5472 msg->
put(maxHydrogenGroupSize);
5473 msg->
put(numMigrationGroups);
5474 msg->
put(maxMigrationGroupSize);
5475 msg->
put(isOccupancyValid);
5476 msg->
put(isBFactorValid);
5479 msg->
put(atomNamePoolSize);
5480 for(
int i=0; i<atomNamePoolSize;i++) {
5487 int numFixedAtomsSet = fixedAtomsSet->size();
5488 msg->
put(numFixedAtoms);
5489 msg->
put(numFixedAtomsSet);
5490 msg->
put(numFixedAtomsSet*
sizeof(AtomSet), (
char *)(fixedAtomsSet->begin()));
5494 int numConstrainedAtomsSet = constrainedAtomsSet->size();
5495 msg->
put(numConstraints);
5496 msg->
put(numConstrainedAtomsSet);
5497 msg->
put(numConstrainedAtomsSet*
sizeof(AtomSet), (
char *)(constrainedAtomsSet->begin()));
5502 msg->
put(numAtoms*
sizeof(
Atom), (
char*)atoms);
5505 msg->
put(numRealBonds);
5510 msg->
put(numBonds*
sizeof(
Bond), (
char*)bonds);
5514 msg->
put(numAngles);
5517 msg->
put(numAngles*
sizeof(
Angle), (
char*)angles);
5521 msg->
put(numDihedrals);
5524 msg->
put(numDihedrals*
sizeof(
Dihedral), (
char*)dihedrals);
5528 msg->
put(num_alch_unpert_Bonds);
5529 msg->
put(num_alch_unpert_Bonds*
sizeof(
Bond), (
char*)alch_unpert_bonds);
5531 msg->
put(num_alch_unpert_Angles);
5532 msg->
put(num_alch_unpert_Angles*
sizeof(
Angle), (
char*)alch_unpert_angles);
5534 msg->
put(num_alch_unpert_Dihedrals);
5535 msg->
put(num_alch_unpert_Dihedrals*
sizeof(
Dihedral), (
char*)alch_unpert_dihedrals);
5539 msg->
put(numImpropers);
5542 msg->
put(numImpropers*
sizeof(
Improper), (
char*)impropers);
5546 msg->
put(numCrossterms);
5549 msg->
put(numCrossterms*
sizeof(
Crossterm), (
char*)crossterms);
5553 msg->
put(numDonors);
5556 msg->
put(numDonors*
sizeof(
Bond), (
char*)donors);
5560 msg->
put(numAcceptors);
5563 msg->
put(numAcceptors*
sizeof(
Bond), (
char*)acceptors);
5567 msg->
put(numExclusions);
5570 msg->
put(numExclusions*
sizeof(
Exclusion), (
char*)exclusions);
5575 msg->
put(numConstraints);
5577 msg->
put(numAtoms, consIndexes);
5581 msg->
put(numConstraints*
sizeof(ConstraintParams), (
char*)consParams);
5590 DebugM(3,
"Sending gridforce info\n" <<
endi);
5591 msg->
put(numGridforceGrids);
5593 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
5594 msg->
put(numGridforces[gridnum]);
5595 msg->
put(numAtoms, gridfrcIndexes[gridnum]);
5596 if (numGridforces[gridnum])
5598 msg->
put(numGridforces[gridnum]*
sizeof(GridforceParams), (
char*)gridfrcParams[gridnum]);
5609 msg->
put(numStirredAtoms);
5611 msg->
put(numAtoms, stirIndexes);
5613 if (numStirredAtoms)
5616 msg->
put(numStirredAtoms*
sizeof(StirParams), (
char*)stirParams);
5623 msg->
put(numMovDrag);
5624 msg->
put(numAtoms, movDragIndexes);
5627 msg->
put(numMovDrag*
sizeof(MovDragParams), (
char*)movDragParams);
5633 msg->
put(numRotDrag);
5634 msg->
put(numAtoms, rotDragIndexes);
5637 msg->
put(numRotDrag*
sizeof(RotDragParams), (
char*)rotDragParams);
5643 msg->
put(numConsTorque);
5644 msg->
put(numAtoms, consTorqueIndexes);
5647 msg->
put(numConsTorque*
sizeof(ConsTorqueParams), (
char*)consTorqueParams);
5653 { msg->
put(numConsForce);
5654 msg->
put(numAtoms, consForceIndexes);
5656 msg->
put(numConsForce*
sizeof(
Vector), (
char*)consForce);
5660 msg->
put(numExPressureAtoms);
5661 msg->
put(numAtoms, exPressureAtomFlags);
5664 #ifndef MEM_OPT_VERSION
5668 msg->
put(numAtoms, langevinParams);
5674 msg->
put(numFixedAtoms);
5675 msg->
put(numAtoms, fixedAtomFlags);
5676 msg->
put(numFixedRigidBonds);
5681 msg->
put(numAtoms, qmAtomGroup);
5682 msg->
put(qmNumQMAtoms);
5683 msg->
put(qmNumQMAtoms, qmAtmChrg);
5684 msg->
put(qmNumQMAtoms, qmAtmIndx);
5686 msg->
put(qmNumBonds);
5687 msg->
put(qmMeNumBonds);
5688 msg->
put(qmMeNumBonds, qmMeMMindx);
5689 msg->
put(qmMeNumBonds, qmMeQMGrp);
5691 msg->
put(qmNumGrps);
5692 msg->
put(qmNumGrps, qmGrpID);
5693 msg->
put(qmNumGrps, qmCustPCSizes);
5694 msg->
put(qmTotCustPCs);
5695 msg->
put(qmTotCustPCs, qmCustomPCIdxs);
5701 msg->
put(numFepInitial);
5702 msg->
put(numFepFinal);
5703 msg->
put(numAtoms*
sizeof(
char), (
char*)fepAtomFlags);
5708 msg->
put(numAtoms*
sizeof(
char), (
char*)ssAtomFlags);
5709 msg->
put(ss_num_vdw_params);
5711 msg->
put(numAtoms*
sizeof(
int), (
char*)ss_index);
5714 #ifdef OPENATOM_VERSION
5716 if (simParams->openatomOn ) {
5717 msg->
put(numFepInitial);
5718 msg->
put(numAtoms*
sizeof(
char), (
char*)fepAtomFlags);
5720 #endif //OPENATOM_VERSION
5723 msg->
put(is_lonepairs_psf);
5724 if (is_lonepairs_psf) {
5725 msg->
put(numLphosts);
5726 msg->
put(numLphosts*
sizeof(
Lphost), (
char*)lphosts);
5728 msg->
put(is_drude_psf);
5731 msg->
put(numTholes);
5732 msg->
put(numTholes*
sizeof(
Thole), (
char*)tholes);
5733 msg->
put(numAnisos);
5734 msg->
put(numAnisos*
sizeof(
Aniso), (
char*)anisos);
5736 msg->
put(numZeroMassAtoms);
5741 msg->
put(numAtoms, (
int*)lcpoParamType);
5746 msg->
put(numLJPair);
5747 msg->
put(numLJPair,indxLJA);
5748 msg->
put(numLJPair,indxLJB);
5749 msg->
put(numLJPair,pairC6);
5750 msg->
put(numLJPair,pairC12);
5751 msg->
put(numLJPair,gromacsPair_type);
5752 msg->
put((numAtoms),pointerToLJBeg);
5753 msg->
put((numAtoms),pointerToLJEnd);
5763 msg->
put((numAtoms),pointerToGaussBeg);
5764 msg->
put((numAtoms),pointerToGaussEnd);
5772 #ifdef MEM_OPT_VERSION
5774 build_excl_check_signatures();
5777 numBonds = numCalcBonds = 0;
5778 numAngles = numCalcAngles = 0;
5779 numDihedrals = numCalcDihedrals = 0;
5780 numImpropers = numCalcImpropers = 0;
5781 numCrossterms = numCalcCrossterms = 0;
5782 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;
5784 numLJPair = numCalcLJPair = 0;
5790 build_lists_by_atom();
5810 #ifdef MEM_OPT_VERSION
5814 msg->
get(massPoolSize);
5815 if(atomMassPool)
delete [] atomMassPool;
5816 atomMassPool =
new Real[massPoolSize];
5817 msg->
get(massPoolSize, atomMassPool);
5819 msg->
get(chargePoolSize);
5820 if(atomChargePool)
delete [] atomChargePool;
5821 atomChargePool =
new Real[chargePoolSize];
5822 msg->
get(chargePoolSize, atomChargePool);
5825 msg->
get(atomSigPoolSize);
5828 for(
int i=0; i<atomSigPoolSize; i++)
5832 msg->
get(exclSigPoolSize);
5833 if(exclSigPool)
delete [] exclSigPool;
5835 for(
int i=0; i<exclSigPoolSize; i++)
5836 exclSigPool[i].unpack(msg);
5838 msg->
get(numHydrogenGroups);
5839 msg->
get(maxHydrogenGroupSize);
5840 msg->
get(numMigrationGroups);
5841 msg->
get(maxMigrationGroupSize);
5842 msg->
get(isOccupancyValid);
5843 msg->
get(isBFactorValid);
5846 msg->
get(atomNamePoolSize);
5848 for(
int i=0; i<atomNamePoolSize;i++) {
5856 int numFixedAtomsSet;
5857 msg->
get(numFixedAtoms);
5858 msg->
get(numFixedAtomsSet);
5859 fixedAtomsSet =
new AtomSetList(numFixedAtomsSet);
5860 msg->
get(numFixedAtomsSet*
sizeof(AtomSet), (
char *)(fixedAtomsSet->begin()));
5864 int numConstrainedAtomsSet;
5865 msg->
get(numConstraints);
5866 msg->
get(numConstrainedAtomsSet);
5867 constrainedAtomsSet =
new AtomSetList(numConstrainedAtomsSet);
5868 msg->
get(numConstrainedAtomsSet*
sizeof(AtomSet), (
char *)(constrainedAtomsSet->begin()));
5873 atoms=
new Atom[numAtoms];
5874 msg->
get(numAtoms*
sizeof(
Atom), (
char*)atoms);
5877 msg->
get(numRealBonds);
5882 bonds=
new Bond[numBonds];
5883 msg->
get(numBonds*
sizeof(
Bond), (
char*)bonds);
5887 msg->
get(numAngles);
5891 angles=
new Angle[numAngles];
5892 msg->
get(numAngles*
sizeof(
Angle), (
char*)angles);
5896 msg->
get(numDihedrals);
5899 delete [] dihedrals;
5900 dihedrals=
new Dihedral[numDihedrals];
5901 msg->
get(numDihedrals*
sizeof(
Dihedral), (
char*)dihedrals);
5905 msg->
get(num_alch_unpert_Bonds);
5906 alch_unpert_bonds=
new Bond[num_alch_unpert_Bonds];
5907 msg->
get(num_alch_unpert_Bonds*
sizeof(
Bond), (
char*)alch_unpert_bonds);
5909 msg->
get(num_alch_unpert_Angles);
5910 alch_unpert_angles=
new Angle[num_alch_unpert_Angles];
5911 msg->
get(num_alch_unpert_Angles*
sizeof(
Angle), (
char*)alch_unpert_angles);
5913 msg->
get(num_alch_unpert_Dihedrals);
5914 alch_unpert_dihedrals=
new Dihedral[num_alch_unpert_Dihedrals];
5915 msg->
get(num_alch_unpert_Dihedrals*
sizeof(
Dihedral), (
char*)alch_unpert_dihedrals);
5919 msg->
get(numImpropers);
5922 delete [] impropers;
5923 impropers=
new Improper[numImpropers];
5924 msg->
get(numImpropers*
sizeof(
Improper), (
char*)impropers);
5928 msg->
get(numCrossterms);
5931 delete [] crossterms;
5932 crossterms=
new Crossterm[numCrossterms];
5933 msg->
get(numCrossterms*
sizeof(
Crossterm), (
char*)crossterms);
5937 msg->
get(numDonors);
5941 donors=
new Bond[numDonors];
5942 msg->
get(numDonors*
sizeof(
Bond), (
char*)donors);
5946 msg->
get(numAcceptors);
5949 delete [] acceptors;
5950 acceptors=
new Bond[numAcceptors];
5951 msg->
get(numAcceptors*
sizeof(
Bond), (
char*)acceptors);
5955 msg->
get(numExclusions);
5959 exclusions=
new Exclusion[numExclusions];
5960 msg->
get(numExclusions*
sizeof(
Exclusion), (
char*)exclusions);
5966 msg->
get(numConstraints);
5968 delete [] consIndexes;
5969 consIndexes =
new int32[numAtoms];
5971 msg->
get(numAtoms, consIndexes);
5975 delete [] consParams;
5976 consParams =
new ConstraintParams[numConstraints];
5978 msg->
get(numConstraints*
sizeof(ConstraintParams), (
char*)consParams);
5986 DebugM(3,
"Receiving gridforce info\n");
5988 msg->
get(numGridforceGrids);
5990 DebugM(3,
"numGridforceGrids = " << numGridforceGrids <<
"\n");
5992 delete [] numGridforces;
5993 numGridforces =
new int[numGridforceGrids];
5995 delete [] gridfrcIndexes;
5996 delete [] gridfrcParams;
5997 delete [] gridfrcGrid;
5998 gridfrcIndexes =
new int32*[numGridforceGrids];
5999 gridfrcParams =
new GridforceParams*[numGridforceGrids];
6002 int grandTotalGrids = 0;
6003 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6004 msg->
get(numGridforces[gridnum]);
6006 gridfrcIndexes[gridnum] =
new int32[numAtoms];
6007 msg->
get(numAtoms, gridfrcIndexes[gridnum]);
6009 if (numGridforces[gridnum])
6011 gridfrcParams[gridnum] =
new GridforceParams[numGridforces[gridnum]];
6012 msg->
get(numGridforces[gridnum]*
sizeof(GridforceParams), (
char*)gridfrcParams[gridnum]);
6017 grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
6025 msg->
get(numStirredAtoms);
6027 delete [] stirIndexes;
6028 stirIndexes =
new int32[numAtoms];
6030 msg->
get(numAtoms, stirIndexes);
6032 if (numStirredAtoms)
6034 delete [] stirParams;
6035 stirParams =
new StirParams[numStirredAtoms];
6037 msg->
get(numStirredAtoms*
sizeof(StirParams), (
char*)stirParams);
6043 msg->
get(numMovDrag);
6044 delete [] movDragIndexes;
6045 movDragIndexes =
new int32[numAtoms];
6046 msg->
get(numAtoms, movDragIndexes);
6049 delete [] movDragParams;
6050 movDragParams =
new MovDragParams[numMovDrag];
6051 msg->
get(numMovDrag*
sizeof(MovDragParams), (
char*)movDragParams);
6057 msg->
get(numRotDrag);
6058 delete [] rotDragIndexes;
6059 rotDragIndexes =
new int32[numAtoms];
6060 msg->
get(numAtoms, rotDragIndexes);
6063 delete [] rotDragParams;
6064 rotDragParams =
new RotDragParams[numRotDrag];
6065 msg->
get(numRotDrag*
sizeof(RotDragParams), (
char*)rotDragParams);
6071 msg->
get(numConsTorque);
6072 delete [] consTorqueIndexes;
6073 consTorqueIndexes =
new int32[numAtoms];
6074 msg->
get(numAtoms, consTorqueIndexes);
6077 delete [] consTorqueParams;
6078 consTorqueParams =
new ConsTorqueParams[numConsTorque];
6079 msg->
get(numConsTorque*
sizeof(ConsTorqueParams), (
char*)consTorqueParams);
6085 { msg->
get(numConsForce);
6086 delete [] consForceIndexes;
6087 consForceIndexes =
new int32[numAtoms];
6088 msg->
get(numAtoms, consForceIndexes);
6090 {
delete [] consForce;
6091 consForce =
new Vector[numConsForce];
6092 msg->
get(numConsForce*
sizeof(
Vector), (
char*)consForce);
6097 exPressureAtomFlags =
new int32[numAtoms];
6098 msg->
get(numExPressureAtoms);
6099 msg->
get(numAtoms, exPressureAtomFlags);
6102 #ifndef MEM_OPT_VERSION
6106 delete [] langevinParams;
6107 langevinParams =
new Real[numAtoms];
6109 msg->
get(numAtoms, langevinParams);
6115 delete [] fixedAtomFlags;
6116 fixedAtomFlags =
new int32[numAtoms];
6118 msg->
get(numFixedAtoms);
6119 msg->
get(numAtoms, fixedAtomFlags);
6120 msg->
get(numFixedRigidBonds);
6125 if( qmAtomGroup != 0)
6126 delete [] qmAtomGroup;
6127 qmAtomGroup =
new Real[numAtoms];
6129 msg->
get(numAtoms, qmAtomGroup);
6131 msg->
get(qmNumQMAtoms);
6134 delete [] qmAtmChrg;
6135 qmAtmChrg =
new Real[qmNumQMAtoms];
6137 msg->
get(qmNumQMAtoms, qmAtmChrg);
6140 delete [] qmAtmIndx;
6141 qmAtmIndx =
new int[qmNumQMAtoms];
6143 msg->
get(qmNumQMAtoms, qmAtmIndx);
6147 msg->
get(qmNumBonds);
6149 msg->
get(qmMeNumBonds);
6151 if( qmMeMMindx != 0)
6152 delete [] qmMeMMindx;
6153 qmMeMMindx =
new int[qmMeNumBonds];
6155 msg->
get(qmMeNumBonds, qmMeMMindx);
6158 delete [] qmMeQMGrp;
6159 qmMeQMGrp =
new Real[qmMeNumBonds];
6161 msg->
get(qmMeNumBonds, qmMeQMGrp);
6165 msg->
get(qmNumGrps);
6169 qmGrpID =
new Real[qmNumGrps];
6170 msg->
get(qmNumGrps, qmGrpID);
6172 if( qmCustPCSizes != 0)
6173 delete [] qmCustPCSizes;
6174 qmCustPCSizes =
new int[qmNumGrps];
6175 msg->
get(qmNumGrps, qmCustPCSizes);
6177 msg->
get(qmTotCustPCs);
6179 if( qmCustomPCIdxs != 0)
6180 delete [] qmCustomPCIdxs;
6181 qmCustomPCIdxs =
new int[qmTotCustPCs];
6182 msg->
get(qmTotCustPCs, qmCustomPCIdxs);
6188 delete [] fepAtomFlags;
6189 fepAtomFlags =
new unsigned char[numAtoms];
6191 msg->
get(numFepInitial);
6192 msg->
get(numFepFinal);
6193 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)fepAtomFlags);
6199 delete [] ssAtomFlags;
6200 delete [] ss_vdw_type;
6202 ssAtomFlags =
new unsigned char[numAtoms];
6204 ss_index =
new int [numAtoms];
6205 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)ssAtomFlags);
6206 msg->
get(ss_num_vdw_params);
6208 msg->
get(numAtoms*
sizeof(
int), (
char*)ss_index);
6211 #ifdef OPENATOM_VERSION
6213 if (simParams->openatomOn) {
6214 delete [] fepAtomFlags;
6215 fepAtomFlags =
new unsigned char[numAtoms];
6217 msg->
get(numFepInitial);
6218 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)fepAtomFlags);
6219 #endif //OPENATOM_VERSION
6222 msg->
get(is_lonepairs_psf);
6223 if (is_lonepairs_psf) {
6224 msg->
get(numLphosts);
6226 lphosts =
new Lphost[numLphosts];
6227 msg->
get(numLphosts*
sizeof(
Lphost), (
char*)lphosts);
6229 msg->
get(is_drude_psf);
6231 delete[] drudeConsts;
6234 msg->
get(numTholes);
6236 tholes =
new Thole[numTholes];
6237 msg->
get(numTholes*
sizeof(
Thole), (
char*)tholes);
6238 msg->
get(numAnisos);
6240 anisos =
new Aniso[numAnisos];
6241 msg->
get(numAnisos*
sizeof(
Aniso), (
char*)anisos);
6243 msg->
get(numZeroMassAtoms);
6248 delete [] lcpoParamType;
6249 lcpoParamType =
new int[numAtoms];
6250 msg->
get(numAtoms, (
int*)lcpoParamType);
6256 msg->
get(numLJPair);
6259 msg->
get(numLJPair,indxLJA);
6262 msg->
get(numLJPair,indxLJB);
6265 msg->
get(numLJPair,pairC6);
6268 msg->
get(numLJPair,pairC12);
6269 delete [] gromacsPair_type;
6271 msg->
get(numLJPair,gromacsPair_type);
6272 delete [] pointerToLJBeg;
6273 pointerToLJBeg =
new int[numAtoms];
6274 msg->
get((numAtoms),pointerToLJBeg);
6275 delete [] pointerToLJEnd;
6276 pointerToLJEnd =
new int[numAtoms];
6277 msg->
get((numAtoms),pointerToLJEnd);
6290 delete [] indxGaussA;
6293 delete [] indxGaussB;
6311 delete [] gRepulsive;
6314 delete [] pointerToGaussBeg;
6315 pointerToGaussBeg =
new int[numAtoms];
6316 msg->
get((numAtoms),pointerToGaussBeg);
6317 delete [] pointerToGaussEnd;
6318 pointerToGaussEnd =
new int[numAtoms];
6319 msg->
get((numAtoms),pointerToGaussEnd);
6327 #ifdef MEM_OPT_VERSION
6329 build_excl_check_signatures();
6332 numBonds = numCalcBonds = 0;
6333 numAngles = numCalcAngles = 0;
6334 numDihedrals = numCalcDihedrals = 0;
6335 numImpropers = numCalcImpropers = 0;
6336 numCrossterms = numCalcCrossterms = 0;
6337 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;
6339 numLJPair = numCalcLJPair = 0;
6345 build_atom_status();
6346 build_lists_by_atom();
6385 DebugM(3,
"Entered build_gridforce_params multi...\n");
6390 numGridforceGrids = 0;
6391 while (mgridParams != NULL) {
6392 numGridforceGrids++;
6393 mgridParams = mgridParams->
next;
6396 DebugM(3,
"numGridforceGrids = " << numGridforceGrids <<
"\n");
6397 gridfrcIndexes =
new int32*[numGridforceGrids];
6398 gridfrcParams =
new GridforceParams*[numGridforceGrids];
6400 numGridforces =
new int[numGridforceGrids];
6402 int grandTotalGrids = 0;
6405 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6406 int current_index=0;
6413 if (mgridParams == NULL) {
6414 NAMD_die(
"Problem with mgridParams!");
6420 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, gridforceFile required.");
6427 if ( (cwd == NULL) || (mgridParams->
gridforceFile[0] ==
'/') )
6433 strcpy(filename, cwd);
6437 kPDB =
new PDB(filename);
6440 NAMD_die(
"Memory allocation failed in Molecule::build_gridforce_params");
6445 NAMD_die(
"Number of atoms in grid force PDB doesn't match coordinate PDB");
6464 else if (strcasecmp(mgridParams->
gridforceCol,
"Y") == 0)
6468 else if (strcasecmp(mgridParams->
gridforceCol,
"Z") == 0)
6472 else if (strcasecmp(mgridParams->
gridforceCol,
"O") == 0)
6476 else if (strcasecmp(mgridParams->
gridforceCol,
"B") == 0)
6482 NAMD_die(
"gridforcecol must have value of X, Y, Z, O, or B");
6515 NAMD_die(
"gridforcechargecol must have value of X, Y, Z, O, or B");
6520 NAMD_die(
"gridforcecol and gridforcechargecol cannot have same value");
6527 gridfrcIndexes[gridnum] =
new int32[numAtoms];
6529 if (gridfrcIndexes[gridnum] == NULL)
6531 NAMD_die(
"memory allocation failed in Molecule::build_gridforce_params()");
6535 for (i=0; i<numAtoms; i++)
6541 kval = (kPDB->
atom(i))->xcoor();
6544 kval = (kPDB->
atom(i))->ycoor();
6547 kval = (kPDB->
atom(i))->zcoor();
6550 kval = (kPDB->
atom(i))->occupancy();
6553 kval = (kPDB->
atom(i))->temperaturefactor();
6560 gridfrcIndexes[gridnum][i] = current_index;
6566 gridfrcIndexes[gridnum][i] = -1;
6570 if (current_index == 0)
6573 iout <<
iWARN <<
"NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" <<
endi;
6578 gridfrcParams[gridnum] =
new GridforceParams[current_index];
6579 if (gridfrcParams[gridnum] == NULL)
6581 NAMD_die(
"memory allocation failed in Molecule::build_gridforce_params");
6585 numGridforces[gridnum] = current_index;
6589 for (i=0; i<numAtoms; i++)
6591 if (gridfrcIndexes[gridnum][i] != -1)
6597 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->xcoor();
6600 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->ycoor();
6603 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->zcoor();
6606 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->occupancy();
6609 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->temperaturefactor();
6617 #ifdef MEM_OPT_VERSION
6618 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
6620 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
6624 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->xcoor();
6627 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->ycoor();
6630 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->zcoor();
6633 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->occupancy();
6636 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->temperaturefactor();
6657 strcpy(potfilename, cwd);
6664 DebugM(3,
"allocating GridforceGrid(" << gridnum <<
")\n");
6667 grandTotalGrids += gridfrcGrid[gridnum]->get_total_grids();
6668 DebugM(4,
"grandTotalGrids = " << grandTotalGrids <<
"\n" <<
endi);
6671 mgridParams = mgridParams->
next;
6677 #endif // MOLECULE2_C undefined = first object file
6678 #ifdef MOLECULE2_C // second object file
6710 int current_index=0;
6720 if (consref == NULL)
6722 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, consref required.");
6723 refPDB = initial_pdb;
6727 if (consref->
next != NULL)
6729 NAMD_die(
"Multiple definitions of constraint reference file in configruation file");
6732 if ( (cwd == NULL) || (consref->
data[0] ==
'/') )
6734 strcpy(filename, consref->
data);
6738 strcpy(filename, cwd);
6739 strcat(filename, consref->
data);
6742 refPDB =
new PDB(filename);
6743 if ( refPDB == NULL )
6745 NAMD_die(
"Memory allocation failed in Molecule::build_constraint_params");
6750 NAMD_die(
"Number of atoms in constraint reference PDB doesn't match coordinate PDB");
6757 if (conskfile == NULL)
6759 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, conskfile required.");
6764 if (conskfile->
next != NULL)
6766 NAMD_die(
"Multiple definitions of constraint constant file in configuration file");
6769 if ( (consref != NULL) && (strcasecmp(consref->
data, conskfile->
data) == 0) )
6776 if ( (cwd == NULL) || (conskfile->
data[0] ==
'/') )
6778 strcpy(filename, conskfile->
data);
6782 strcpy(filename, cwd);
6783 strcat(filename, conskfile->
data);
6786 kPDB =
new PDB(filename);
6789 NAMD_die(
"Memory allocation failed in Molecule::build_constraint_params");
6794 NAMD_die(
"Number of atoms in constraint constant PDB doesn't match coordinate PDB");
6804 if (conskcol == NULL)
6810 if (conskcol->
next != NULL)
6812 NAMD_die(
"Multiple definitions of harmonic constraint column in config file");
6815 if (strcasecmp(conskcol->
data,
"X") == 0)
6819 else if (strcasecmp(conskcol->
data,
"Y") == 0)
6823 else if (strcasecmp(conskcol->
data,
"Z") == 0)
6827 else if (strcasecmp(conskcol->
data,
"O") == 0)
6831 else if (strcasecmp(conskcol->
data,
"B") == 0)
6837 NAMD_die(
"conskcol must have value of X, Y, Z, O, or B");
6844 consIndexes =
new int32[numAtoms];
6846 if (consIndexes == NULL)
6848 NAMD_die(
"memory allocation failed in Molecule::build_constraint_params()");
6852 for (i=0; i<numAtoms; i++)
6858 kval = (kPDB->
atom(i))->xcoor();
6861 kval = (kPDB->
atom(i))->ycoor();
6864 kval = (kPDB->
atom(i))->zcoor();
6867 kval = (kPDB->
atom(i))->occupancy();
6870 kval = (kPDB->
atom(i))->temperaturefactor();
6877 consIndexes[i] = current_index;
6883 consIndexes[i] = -1;
6887 if (current_index == 0)
6890 iout <<
iWARN <<
"NO CONSTRAINED ATOMS WERE FOUND, BUT CONSTRAINTS ARE ON . . .\n" <<
endi;
6895 consParams =
new ConstraintParams[current_index];
6897 if (consParams == NULL)
6899 NAMD_die(
"memory allocation failed in Molecule::build_constraint_params");
6903 numConstraints = current_index;
6907 for (i=0; i<numAtoms; i++)
6909 if (consIndexes[i] != -1)
6915 consParams[consIndexes[i]].k = (kPDB->
atom(i))->xcoor();
6918 consParams[consIndexes[i]].k = (kPDB->
atom(i))->ycoor();
6921 consParams[consIndexes[i]].k = (kPDB->
atom(i))->zcoor();
6924 consParams[consIndexes[i]].k = (kPDB->
atom(i))->occupancy();
6927 consParams[consIndexes[i]].k = (kPDB->
atom(i))->temperaturefactor();
6932 consParams[consIndexes[i]].refPos.x = (refPDB->
atom(i))->xcoor();
6933 consParams[consIndexes[i]].refPos.y = (refPDB->
atom(i))->ycoor();
6934 consParams[consIndexes[i]].refPos.z = (refPDB->
atom(i))->zcoor();
6939 if (consref != NULL)
6944 if ((conskfile != NULL) &&
6945 !((consref != NULL) &&
6946 (strcasecmp(consref->
data, conskfile->
data) == 0)
6985 int current_index=0;
6994 if (movDragFile == NULL) {
6995 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, movDragFile required.");
7000 if (movDragFile->
next != NULL) {
7001 NAMD_die(
"Multiple definitions of moving drag tag file in configuration file");
7004 if ( (cwd == NULL) || (movDragFile->
data[0] ==
'/') ) {
7005 strcpy(mainfilename, movDragFile->
data);
7007 strcpy(mainfilename, cwd);
7008 strcat(mainfilename, movDragFile->
data);
7011 tPDB =
new PDB(mainfilename);
7012 if ( tPDB == NULL ) {
7013 NAMD_die(
"Memory allocation failed in Molecule::build_movdrag_params");
7017 NAMD_die(
"Number of atoms in moving drag tag PDB doesn't match coordinate PDB");
7025 if (movDragVelFile == NULL) {
7026 if (movDragFile == NULL) {
7027 NAMD_die(
"Moving drag velocity file can not be same as coordinate PDB file");
7029 if (movDragVelFile->
next != NULL) {
7030 NAMD_die(
"Multiple definitions of moving drag velocity file in configuration file");
7037 if ( (cwd == NULL) || (movDragVelFile->
data[0] ==
'/') ) {
7038 strcpy(velfilename, movDragVelFile->
data);
7040 strcpy(velfilename, cwd);
7041 strcat(velfilename, movDragVelFile->
data);
7044 vPDB =
new PDB(velfilename);
7045 if ( vPDB == NULL ) {
7046 NAMD_die(
"Memory allocation failed in Molecule::build_movdrag_params");
7050 NAMD_die(
"Number of atoms in moving drag velocity PDB doesn't match coordinate PDB");
7062 if (movDragCol == NULL) {
7065 if (movDragCol->
next != NULL) {
7066 NAMD_die(
"Multiple definitions of drag column in config file");
7069 if (movDragFile == NULL
7070 && strcasecmp(movDragCol->
data,
"B")
7071 && strcasecmp(movDragCol->
data,
"O")) {
7072 NAMD_die(
"Can not read moving drag tags from X, Y, or Z column of the coordinate or velocity file");
7074 if (!strcasecmp(movDragCol->
data,
"X")) {
7076 }
else if (!strcasecmp(movDragCol->
data,
"Y")) {
7078 }
else if (!strcasecmp(movDragCol->
data,
"Z")) {
7080 }
else if (!strcasecmp(movDragCol->
data,
"O")) {
7082 }
else if (!strcasecmp(movDragCol->
data,
"B")) {
7086 NAMD_die(
"movDragCol must have value of X, Y, Z, O, or B");
7093 movDragIndexes =
new int32[numAtoms];
7094 if (movDragIndexes == NULL) {
7095 NAMD_die(
"memory allocation failed in Molecule::build_movdrag_params()");
7099 for (i=0; i<numAtoms; i++) {
7102 dtval = (tPDB->
atom(i))->xcoor();
7105 dtval = (tPDB->
atom(i))->ycoor();
7108 dtval = (tPDB->
atom(i))->zcoor();
7111 dtval = (tPDB->
atom(i))->occupancy();
7114 dtval = (tPDB->
atom(i))->temperaturefactor();
7120 movDragIndexes[i] = current_index;
7124 movDragIndexes[i] = -1;
7128 if (current_index == 0) {
7130 iout <<
iWARN <<
"NO DRAGGED ATOMS WERE FOUND, BUT MOVING DRAG IS ON . . . " <<
endi;
7133 movDragParams =
new MovDragParams[current_index];
7134 if (movDragParams == NULL) {
7135 NAMD_die(
"memory allocation failed in Molecule::build_movdrag_params");
7139 numMovDrag = current_index;
7143 for (i=0; i<numAtoms; i++) {
7144 if (movDragIndexes[i] != -1) {
7145 movDragParams[movDragIndexes[i]].v[0] = (vPDB->
atom(i))->xcoor();
7146 movDragParams[movDragIndexes[i]].v[1] = (vPDB->
atom(i))->ycoor();
7147 movDragParams[movDragIndexes[i]].v[2] = (vPDB->
atom(i))->zcoor();
7151 if (movDragFile != NULL)
delete tPDB;
7152 if (movDragVelFile != NULL)
delete vPDB;
7189 PDB *tPDB, *aPDB, *pPDB, *vPDB;
7191 int current_index=0;
7204 if (rotDragFile == NULL) {
7205 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, rotDragFile required.");
7210 if (rotDragFile->
next != NULL) {
7211 NAMD_die(
"Multiple definitions of rotating drag tag file in configuration file");
7214 if ( (cwd == NULL) || (rotDragFile->
data[0] ==
'/') ) {
7215 strcpy(mainfilename, rotDragFile->
data);
7217 strcpy(mainfilename, cwd);
7218 strcat(mainfilename, rotDragFile->
data);
7221 tPDB =
new PDB(mainfilename);
7222 if ( tPDB == NULL ) {
7223 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7227 NAMD_die(
"Number of atoms in rotating drag tag PDB doesn't match coordinate PDB");
7235 if (rotDragAxisFile == NULL) {
7236 if (rotDragFile == NULL) {
7237 NAMD_die(
"Rotating drag axis file can not be same as coordinate PDB file");
7239 if (rotDragAxisFile->
next != NULL) {
7240 NAMD_die(
"Multiple definitions of rotating drag axis file in configuration file");
7242 if (rotDragPivotFile == NULL) {
7243 NAMD_die(
"Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
7250 if ( (cwd == NULL) || (rotDragAxisFile->
data[0] ==
'/') ) {
7251 strcpy(axisfilename, rotDragAxisFile->
data);
7253 strcpy(axisfilename, cwd);
7254 strcat(axisfilename, rotDragAxisFile->
data);
7257 aPDB =
new PDB(axisfilename);
7258 if ( aPDB == NULL ) {
7259 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7263 NAMD_die(
"Number of atoms in rotating drag axis PDB doesn't match coordinate PDB");
7271 if (rotDragPivotFile == NULL) {
7272 if (rotDragFile == NULL) {
7273 NAMD_die(
"Rotating drag pivot point file can not be same as coordinate PDB file");
7275 if (rotDragPivotFile->
next != NULL) {
7276 NAMD_die(
"Multiple definitions of rotating drag pivot point file in configuration file");
7278 if (rotDragAxisFile == NULL) {
7279 NAMD_die(
"Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
7286 if ( (cwd == NULL) || (rotDragPivotFile->
data[0] ==
'/') ) {
7287 strcpy(pivotfilename, rotDragPivotFile->
data);
7289 strcpy(pivotfilename, cwd);
7290 strcat(pivotfilename, rotDragPivotFile->
data);
7293 pPDB =
new PDB(pivotfilename);
7294 if ( pPDB == NULL ) {
7295 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7299 NAMD_die(
"Number of atoms in rotating drag pivot point PDB doesn't match coordinate PDB");
7308 if (rotDragVelFile == NULL) {
7311 if (rotDragVelFile->
next != NULL) {
7312 NAMD_die(
"Multiple definitions of rotating drag velocity file in configuration file");
7315 if ( (cwd == NULL) || (rotDragVelFile->
data[0] ==
'/') ) {
7316 strcpy(velfilename, rotDragVelFile->
data);
7318 strcpy(velfilename, cwd);
7319 strcat(velfilename, rotDragVelFile->
data);
7322 vPDB =
new PDB(velfilename);
7323 if ( vPDB == NULL ) {
7324 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7328 NAMD_die(
"Number of atoms in rotating drag velocity PDB doesn't match coordinate PDB");
7339 if (rotDragCol == NULL) {
7342 if (rotDragCol->
next != NULL) {
7343 NAMD_die(
"Multiple definitions of drag tag column in config file");
7346 if ( rotDragFile == NULL
7347 && (!strcasecmp(rotDragCol->
data,
"X")
7348 || !strcasecmp(rotDragCol->
data,
"Y")
7349 || !strcasecmp(rotDragCol->
data,
"Z"))) {
7350 NAMD_die(
"Can not read rotating drag tags from X, Y, or Z column of the PDB coordinate file");
7352 if (!strcasecmp(rotDragCol->
data,
"X")) {
7354 }
else if (!strcasecmp(rotDragCol->
data,
"Y")) {
7356 }
else if (!strcasecmp(rotDragCol->
data,
"Z")) {
7358 }
else if (!strcasecmp(rotDragCol->
data,
"O")) {
7360 }
else if (!strcasecmp(rotDragCol->
data,
"B")) {
7364 NAMD_die(
"rotDragCol must have value of X, Y, Z, O, or B");
7376 if (rotDragVelCol == NULL) {
7379 if (rotDragVelCol->
next != NULL) {
7380 NAMD_die(
"Multiple definitions of drag angular velocity column in config file");
7383 if (rotDragVelFile == NULL
7384 && rotDragFile == NULL
7385 && strcasecmp(rotDragCol->
data,
"B")
7386 && strcasecmp(rotDragCol->
data,
"O")) {
7387 NAMD_die(
"Can not read rotating drag angular velocities from X, Y, or Z column of the PDB coordinate file");
7389 if (!strcasecmp(rotDragVelCol->
data,
"X")) {
7391 }
else if (!strcasecmp(rotDragVelCol->
data,
"Y")) {
7393 }
else if (!strcasecmp(rotDragVelCol->
data,
"Z")) {
7395 }
else if (!strcasecmp(rotDragVelCol->
data,
"O")) {
7397 }
else if (!strcasecmp(rotDragVelCol->
data,
"B")) {
7401 NAMD_die(
"rotDragVelCol must have value of X, Y, Z, O, or B");
7408 rotDragIndexes =
new int32[numAtoms];
7409 if (rotDragIndexes == NULL) {
7410 NAMD_die(
"memory allocation failed in Molecule::build_rotdrag_params()");
7414 for (i=0; i<numAtoms; i++) {
7417 dtval = (tPDB->
atom(i))->xcoor();
7420 dtval = (tPDB->
atom(i))->ycoor();
7423 dtval = (tPDB->
atom(i))->zcoor();
7426 dtval = (tPDB->
atom(i))->occupancy();
7429 dtval = (tPDB->
atom(i))->temperaturefactor();
7435 rotDragIndexes[i] = current_index;
7439 rotDragIndexes[i] = -1;
7443 if (current_index == 0) {
7444 iout <<
iWARN <<
"NO DRAGGED ATOMS WERE FOUND, BUT ROTATING DRAG IS ON . . . " <<
endi;
7446 rotDragParams =
new RotDragParams[current_index];
7447 if (rotDragParams == NULL) {
7448 NAMD_die(
"memory allocation failed in Molecule::build_rotdrag_params");
7452 numRotDrag = current_index;
7456 for (i=0; i<numAtoms; i++) {
7457 if (rotDragIndexes[i] != -1) {
7458 rotDragParams[rotDragIndexes[i]].a[0] = (aPDB->
atom(i))->xcoor();
7459 rotDragParams[rotDragIndexes[i]].a[1] = (aPDB->
atom(i))->ycoor();
7460 rotDragParams[rotDragIndexes[i]].a[2] = (aPDB->
atom(i))->zcoor();
7461 rotDragParams[rotDragIndexes[i]].p[0] = (pPDB->
atom(i))->xcoor();
7462 rotDragParams[rotDragIndexes[i]].p[1] = (pPDB->
atom(i))->ycoor();
7463 rotDragParams[rotDragIndexes[i]].p[2] = (pPDB->
atom(i))->zcoor();
7466 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->xcoor();
7469 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->ycoor();
7472 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->zcoor();
7475 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->occupancy();
7478 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->temperaturefactor();
7484 if (rotDragFile != NULL)
delete tPDB;
7485 if (rotDragAxisFile != NULL)
delete aPDB;
7486 if (rotDragPivotFile != NULL)
delete pPDB;
7487 if (rotDragVelFile != NULL)
delete vPDB;
7524 PDB *tPDB, *aPDB, *pPDB, *vPDB;
7526 int current_index=0;
7539 if (consTorqueFile == NULL) {
7540 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, consTorqueFile required.");
7545 if (consTorqueFile->
next != NULL) {
7546 NAMD_die(
"Multiple definitions of \"constant\" torque tag file in configuration file");
7549 if ( (cwd == NULL) || (consTorqueFile->
data[0] ==
'/') ) {
7550 strcpy(mainfilename, consTorqueFile->
data);
7552 strcpy(mainfilename, cwd);
7553 strcat(mainfilename, consTorqueFile->
data);
7556 tPDB =
new PDB(mainfilename);
7557 if ( tPDB == NULL ) {
7558 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7562 NAMD_die(
"Number of atoms in \"constant\" torque tag PDB doesn't match coordinate PDB");
7570 if (consTorqueAxisFile == NULL) {
7571 if (consTorqueFile == NULL) {
7572 NAMD_die(
"\"Constant\" torque axis file can not be same as coordinate PDB file");
7574 if (consTorqueAxisFile->
next != NULL) {
7575 NAMD_die(
"Multiple definitions of \"constant\" torque axis file in configuration file");
7577 if (consTorquePivotFile == NULL) {
7578 NAMD_die(
"Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
7585 if ( (cwd == NULL) || (consTorqueAxisFile->
data[0] ==
'/') ) {
7586 strcpy(axisfilename, consTorqueAxisFile->
data);
7588 strcpy(axisfilename, cwd);
7589 strcat(axisfilename, consTorqueAxisFile->
data);
7592 aPDB =
new PDB(axisfilename);
7593 if ( aPDB == NULL ) {
7594 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7598 NAMD_die(
"Number of atoms in \"constant\" torque axis PDB doesn't match coordinate PDB");
7606 if (consTorquePivotFile == NULL) {
7607 if (consTorqueFile == NULL) {
7608 NAMD_die(
"\"Constant\" torque pivot point file can not be same as coordinate PDB file");
7610 if (consTorquePivotFile->
next != NULL) {
7611 NAMD_die(
"Multiple definitions of \"constant\" torque pivot point file in configuration file");
7613 if (consTorqueAxisFile == NULL) {
7614 NAMD_die(
"Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
7621 if ( (cwd == NULL) || (consTorquePivotFile->
data[0] ==
'/') ) {
7622 strcpy(pivotfilename, consTorquePivotFile->
data);
7624 strcpy(pivotfilename, cwd);
7625 strcat(pivotfilename, consTorquePivotFile->
data);
7628 pPDB =
new PDB(pivotfilename);
7629 if ( pPDB == NULL ) {
7630 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7634 NAMD_die(
"Number of atoms in \"constant\" torque pivot point PDB doesn't match coordinate PDB");
7643 if (consTorqueValFile == NULL) {
7646 if (consTorqueValFile->
next != NULL) {
7647 NAMD_die(
"Multiple definitions of \"constant\" torque velocity file in configuration file");
7650 if ( (cwd == NULL) || (consTorqueValFile->
data[0] ==
'/') ) {
7651 strcpy(velfilename, consTorqueValFile->
data);
7653 strcpy(velfilename, cwd);
7654 strcat(velfilename, consTorqueValFile->
data);
7657 vPDB =
new PDB(velfilename);
7658 if ( vPDB == NULL ) {
7659 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7663 NAMD_die(
"Number of atoms in \"constant\" torque velocity PDB doesn't match coordinate PDB");
7674 if (consTorqueCol == NULL) {
7677 if (consTorqueCol->
next != NULL) {
7678 NAMD_die(
"Multiple definitions of torque tag column in config file");
7681 if ( consTorqueFile == NULL
7682 && (!strcasecmp(consTorqueCol->
data,
"X")
7683 || !strcasecmp(consTorqueCol->
data,
"Y")
7684 || !strcasecmp(consTorqueCol->
data,
"Z"))) {
7685 NAMD_die(
"Can not read \"constant\" torque tags from X, Y, or Z column of the PDB coordinate file");
7687 if (!strcasecmp(consTorqueCol->
data,
"X")) {
7689 }
else if (!strcasecmp(consTorqueCol->
data,
"Y")) {
7691 }
else if (!strcasecmp(consTorqueCol->
data,
"Z")) {
7693 }
else if (!strcasecmp(consTorqueCol->
data,
"O")) {
7695 }
else if (!strcasecmp(consTorqueCol->
data,
"B")) {
7699 NAMD_die(
"consTorqueCol must have value of X, Y, Z, O, or B");
7711 if (consTorqueValCol == NULL) {
7714 if (consTorqueValCol->
next != NULL) {
7715 NAMD_die(
"Multiple definitions of torque value column in config file");
7718 if (consTorqueValFile == NULL
7719 && consTorqueFile == NULL
7720 && strcasecmp(consTorqueCol->
data,
"B")
7721 && strcasecmp(consTorqueCol->
data,
"O")) {
7722 NAMD_die(
"Can not read \"constant\" torque values from X, Y, or Z column of the PDB coordinate file");
7724 if (!strcasecmp(consTorqueValCol->
data,
"X")) {
7726 }
else if (!strcasecmp(consTorqueValCol->
data,
"Y")) {
7728 }
else if (!strcasecmp(consTorqueValCol->
data,
"Z")) {
7730 }
else if (!strcasecmp(consTorqueValCol->
data,
"O")) {
7732 }
else if (!strcasecmp(consTorqueValCol->
data,
"B")) {
7736 NAMD_die(
"consTorqueValCol must have value of X, Y, Z, O, or B");
7743 consTorqueIndexes =
new int32[numAtoms];
7744 if (consTorqueIndexes == NULL) {
7745 NAMD_die(
"memory allocation failed in Molecule::build_constorque_params()");
7749 for (i=0; i<numAtoms; i++) {
7752 dtval = (tPDB->
atom(i))->xcoor();
7755 dtval = (tPDB->
atom(i))->ycoor();
7758 dtval = (tPDB->
atom(i))->zcoor();
7761 dtval = (tPDB->
atom(i))->occupancy();
7764 dtval = (tPDB->
atom(i))->temperaturefactor();
7770 consTorqueIndexes[i] = current_index;
7774 consTorqueIndexes[i] = -1;
7778 if (current_index == 0) {
7779 iout <<
iWARN <<
"NO TORQUED ATOMS WERE FOUND, BUT \"CONSTANT\" TORQUE IS ON . . . " <<
endi;
7781 consTorqueParams =
new ConsTorqueParams[current_index];
7782 if (consTorqueParams == NULL) {
7783 NAMD_die(
"memory allocation failed in Molecule::build_constorque_params");
7787 numConsTorque = current_index;
7791 for (i=0; i<numAtoms; i++) {
7792 if (consTorqueIndexes[i] != -1) {
7793 consTorqueParams[consTorqueIndexes[i]].a[0] = (aPDB->
atom(i))->xcoor();
7794 consTorqueParams[consTorqueIndexes[i]].a[1] = (aPDB->
atom(i))->ycoor();
7795 consTorqueParams[consTorqueIndexes[i]].a[2] = (aPDB->
atom(i))->zcoor();
7796 consTorqueParams[consTorqueIndexes[i]].p[0] = (pPDB->
atom(i))->xcoor();
7797 consTorqueParams[consTorqueIndexes[i]].p[1] = (pPDB->
atom(i))->ycoor();
7798 consTorqueParams[consTorqueIndexes[i]].p[2] = (pPDB->
atom(i))->zcoor();
7801 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->xcoor();
7804 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->ycoor();
7807 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->zcoor();
7810 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->occupancy();
7813 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->temperaturefactor();
7819 if (consTorqueFile != NULL)
delete tPDB;
7820 if (consTorqueAxisFile != NULL)
delete aPDB;
7821 if (consTorquePivotFile != NULL)
delete pPDB;
7822 if (consTorqueValFile != NULL)
delete vPDB;
7848 iout <<
iWARN <<
"NO CONSTANT FORCES SPECIFIED, BUT CONSTANT FORCE IS ON . . .\n" <<
endi;
7849 consForceIndexes =
new int32[numAtoms];
7850 for (i=0; i<numAtoms; i++) consForceIndexes[i] = -1;
7854 if ((forcePDB=
new PDB(filename)) == NULL)
7855 NAMD_die(
"Memory allocation failed in Molecule::build_constant_forces");
7857 NAMD_die(
"Number of atoms in constant force PDB doesn't match coordinate PDB");
7862 consForceIndexes =
new int32[numAtoms];
7863 if (consForceIndexes == NULL)
7864 NAMD_die(
"memory allocation failed in Molecule::build_constant_forces()");
7868 for (i=0; i<numAtoms; i++)
7872 consForceIndexes[i] = -1;
7875 consForceIndexes[i] = numConsForce++;
7877 if (numConsForce == 0)
7879 iout <<
iWARN <<
"NO NON-ZERO FORCES WERE FOUND, BUT CONSTANT FORCE IS ON . . .\n" <<
endi;
7882 consForce =
new Vector[numConsForce];
7883 if (consForce == NULL)
7884 NAMD_die(
"memory allocation failed in Molecule::build_constant_forces");
7886 for (i=0; i<numAtoms; i++)
7887 if ((index=consForceIndexes[i]) != -1)
7904 langevinParams =
new Real[numAtoms];
7906 if ( (langevinParams == NULL) )
7908 NAMD_die(
"memory allocation failed in Molecule::build_langevin_params()");
7912 for (
int i=0; i<numAtoms; i++)
7916 if ( (! doHydrogen) && is_hydrogen(i) ) bval = 0;
7917 else if ( is_lp(i) ) bval = 0;
7918 else if ( is_drude(i) ) bval = drudeCoupling;
7921 langevinParams[i] = bval;
7958 if (langfile == NULL)
7960 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, langevinFile required.");
7965 if (langfile->
next != NULL)
7967 NAMD_die(
"Multiple definitions of langvein PDB file in configuration file");
7970 if ( (cwd == NULL) || (langfile->
data[0] ==
'/') )
7972 strcpy(filename, langfile->
data);
7976 strcpy(filename, cwd);
7977 strcat(filename, langfile->
data);
7980 bPDB =
new PDB(filename);
7983 NAMD_die(
"Memory allocation failed in Molecule::build_langevin_params");
7988 NAMD_die(
"Number of atoms in langevin parameter PDB doesn't match coordinate PDB");
7997 if (langcol == NULL)
8003 if (langcol->
next != NULL)
8005 NAMD_die(
"Multiple definitions of langevin parameter column in config file");
8008 if (strcasecmp(langcol->
data,
"X") == 0)
8012 else if (strcasecmp(langcol->
data,
"Y") == 0)
8016 else if (strcasecmp(langcol->
data,
"Z") == 0)
8020 else if (strcasecmp(langcol->
data,
"O") == 0)
8024 else if (strcasecmp(langcol->
data,
"B") == 0)
8030 NAMD_die(
"langevincol must have value of X, Y, Z, O, or B");
8035 langevinParams =
new Real[numAtoms];
8037 if ( (langevinParams == NULL) )
8039 NAMD_die(
"memory allocation failed in Molecule::build_langevin_params()");
8043 for (i=0; i<numAtoms; i++)
8049 bval = (bPDB->
atom(i))->xcoor();
8052 bval = (bPDB->
atom(i))->ycoor();
8055 bval = (bPDB->
atom(i))->zcoor();
8058 bval = (bPDB->
atom(i))->occupancy();
8061 bval = (bPDB->
atom(i))->temperaturefactor();
8066 langevinParams[i] = bval;
8070 if (langfile != NULL)
8109 if (fixedfile == NULL)
8111 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, fixedAtomsFile required.");
8116 if (fixedfile->
next != NULL)
8118 NAMD_die(
"Multiple definitions of fixed atoms PDB file in configuration file");
8121 if ( (cwd == NULL) || (fixedfile->
data[0] ==
'/') )
8123 strcpy(filename, fixedfile->
data);
8127 strcpy(filename, cwd);
8128 strcat(filename, fixedfile->
data);
8131 bPDB =
new PDB(filename);
8134 NAMD_die(
"Memory allocation failed in Molecule::build_fixed_atoms");
8139 NAMD_die(
"Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
8148 if (fixedcol == NULL)
8154 if (fixedcol->
next != NULL)
8156 NAMD_die(
"Multiple definitions of fixed atoms column in config file");
8159 if (strcasecmp(fixedcol->
data,
"X") == 0)
8163 else if (strcasecmp(fixedcol->
data,
"Y") == 0)
8167 else if (strcasecmp(fixedcol->
data,
"Z") == 0)
8171 else if (strcasecmp(fixedcol->
data,
"O") == 0)
8175 else if (strcasecmp(fixedcol->
data,
"B") == 0)
8181 NAMD_die(
"fixedatomscol must have value of X, Y, Z, O, or B");
8186 fixedAtomFlags =
new int32[numAtoms];
8188 if (fixedAtomFlags == NULL)
8190 NAMD_die(
"memory allocation failed in Molecule::build_fixed_atoms()");
8196 for (i=0; i<numAtoms; i++)
8202 bval = (bPDB->
atom(i))->xcoor();
8205 bval = (bPDB->
atom(i))->ycoor();
8208 bval = (bPDB->
atom(i))->zcoor();
8211 bval = (bPDB->
atom(i))->occupancy();
8214 bval = (bPDB->
atom(i))->temperaturefactor();
8220 fixedAtomFlags[i] = 1;
8224 fixedAtomFlags[i] = 0;
8229 if (fixedfile != NULL)
8236 if ( numRigidBonds ) {
8238 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
8239 int parentIsFixed = 0;
8240 for( ; h_i != h_e; ++h_i ) {
8242 parentIsFixed = fixedAtomFlags[h_i->atomID];
8243 if ( (rigidBondLengths[h_i->atomID] > 0.)
8244 && fixedAtomFlags[h_i[1].atomID]
8245 && fixedAtomFlags[h_i[2].atomID] ) {
8246 ++numFixedRigidBonds;
8249 if ( (rigidBondLengths[h_i->atomID] > 0.)
8250 && fixedAtomFlags[h_i->atomID]
8251 && parentIsFixed ) {
8252 ++numFixedRigidBonds;
8262 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
8264 for( ; h_i != h_e; ++h_i ) {
8266 if ( allFixed ) ++numFixedGroups;
8269 allFixed = allFixed && fixedAtomFlags[h_i->atomID];
8271 if ( allFixed ) ++numFixedGroups;
8308 int current_index=0;
8316 if (stirredfile == NULL)
8318 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, stirFilename required.");
8321 iout <<
iWARN <<
"STIRRING USING INITIAL POSITION FILE FOR REFERENCE POSITIONS" <<
endi;
8325 if (stirredfile->
next != NULL)
8327 NAMD_die(
"Multiple definitions of stirred atoms PDB file in configuration file");
8330 if ( (cwd == NULL) || (stirredfile->
data[0] ==
'/') )
8332 strcpy(filename, stirredfile->
data);
8336 strcpy(filename, cwd);
8337 strcat(filename, stirredfile->
data);
8341 sPDB =
new PDB(filename);
8345 NAMD_die(
"Memory allocation failed in Molecule::build_stirred_atoms");
8351 NAMD_die(
"Number of atoms in stirred atoms PDB doesn't match coordinate PDB");
8363 if (stirredcol == NULL)
8369 if (stirredcol->
next != NULL)
8371 NAMD_die(
"Multiple definitions of stirred atoms column in config file");
8374 if (strcasecmp(stirredcol->
data,
"O") == 0)
8378 else if (strcasecmp(stirredcol->
data,
"B") == 0)
8384 NAMD_die(
"stirredAtomsCol must have value of O or B");
8391 stirIndexes =
new int32[numAtoms];
8393 if (stirIndexes == NULL)
8395 NAMD_die(
"memory allocation failed in Molecule::build_stirred_params()");
8400 for (i=0; i<numAtoms; i++)
8410 bval = (sPDB->
atom(i))->occupancy();
8413 bval = (sPDB->
atom(i))->temperaturefactor();
8422 stirIndexes[i] = current_index;
8428 stirIndexes[i] = -1;
8436 if (current_index == 0)
8439 iout <<
iWARN <<
"NO STIRRED ATOMS WERE FOUND, BUT STIRRING TORQUES ARE ON . . .\n" <<
endi;
8444 stirParams =
new StirParams[current_index];
8446 if (stirParams == NULL)
8448 NAMD_die(
"memory allocation failed in Molecule::build_stir_params");
8452 numStirredAtoms = current_index;
8456 for (i=0; i<numAtoms; i++)
8458 if (stirIndexes[i] != -1)
8462 stirParams[stirIndexes[i]].refPos.x = (sPDB->
atom(i))->xcoor();
8463 stirParams[stirIndexes[i]].refPos.y = (sPDB->
atom(i))->ycoor();
8464 stirParams[stirIndexes[i]].refPos.z = (sPDB->
atom(i))->zcoor();
8469 if (stirredfile != NULL)
8485 int a1,a2,a3,a4;
float k, ref, upper;
8487 #ifndef MEM_OPT_VERSION
8500 NAMD_die(
"NO EXTRA BONDS FILES SPECIFIED");
8503 for ( ; file; file = file->
next ) {
8504 FILE *f = fopen(file->
data,
"r");
8506 sprintf(err_msg,
"UNABLE TO OPEN EXTRA BONDS FILE %s", file->
data);
8518 if (ret_code!=0)
break;
8521 sscanf(buffer,
"%s",type);
8523 #define CHECKATOMID(ATOMID) if ( ATOMID < 0 || ATOMID >= numAtoms ) badatom = 1;
8527 if ( ! strncasecmp(type,
"bond",4) ) {
8528 if ( sscanf(buffer,
"%s %d %d %f %f %s",
8529 type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1;
8535 #ifndef MEM_OPT_VERSION
8538 tmp.atom1 = a1; tmp.atom2 = a2;
8543 tmpv.
k = k; tmpv.
x0 = ref;
8544 bond_params.
add(tmpv);
8545 }
else if ( ! strncasecmp(type,
"wall",4) ) {
8548 if ( sscanf(buffer,
"%s %d %d %f %f %f %s",
8549 type, &a1, &a2, &k, &ref, &upper, err_msg) != 6 ) badline = 1;
8550 else if (upper < ref) badline = 1;
8556 #ifndef MEM_OPT_VERSION
8559 tmp.atom1 = a1; tmp.atom2 = a2;
8564 tmpv.
k = k; tmpv.
x0 = ref; tmpv.
x1 = upper;
8565 bond_params.
add(tmpv);
8566 }
else if ( ! strncasecmp(type,
"angle",4) ) {
8567 if ( sscanf(buffer,
"%s %d %d %d %f %f %s",
8568 type, &a1, &a2, &a3, &k, &ref, err_msg) != 6 ) badline = 1;
8574 #ifndef MEM_OPT_VERSION
8576 tmp.
atom1 = a1; tmp.atom2 = a2; tmp.atom3 = a3;
8582 tmpv.
k = k; tmpv.
theta0 = ref / 180. *
PI;
8584 tmpv.
normal = anglesNormal;
8585 angle_params.
add(tmpv);
8587 }
else if ( ! strncasecmp(type,
"dihedral",4) ) {
8589 int ret = 1 + sscanf(buffer,
"%s %d %d %d %d %f %f %s",
8590 type, &a1, &a2, &a3, &a4, &k, &ref, err_msg);
8592 ret = sscanf(buffer,
"%s %d %d %d %d %f %d %f %s",
8593 type, &a1, &a2, &a3, &a4, &k, &n, &ref, err_msg);
8595 if ( ret != 8 ) badline = 1;
8602 #ifndef MEM_OPT_VERSION
8604 tmp.
atom1 = a1; tmp.atom2 = a2; tmp.atom3 = a3; tmp.atom4 = a4;
8612 dihedral_params.
add(tmpv);
8613 }
else if ( ! strncasecmp(type,
"improper",4) ) {
8615 int ret = 1 + sscanf(buffer,
"%s %d %d %d %d %f %f %s",
8616 type, &a1, &a2, &a3, &a4, &k, &ref, err_msg);
8618 ret = sscanf(buffer,
"%s %d %d %d %d %f %d %f %s",
8619 type, &a1, &a2, &a3, &a4, &k, &n, &ref, err_msg);
8621 if ( ret != 8 ) badline = 1;
8628 #ifndef MEM_OPT_VERSION
8630 tmp.
atom1 = a1; tmp.atom2 = a2; tmp.atom3 = a3; tmp.atom4 = a4;
8638 improper_params.
add(tmpv);
8639 }
else if ( ! strncasecmp(type,
"#",1) ) {
8646 sprintf(err_msg,
"BAD LINE IN EXTRA BONDS FILE %s: %s",
8647 file->
data, buffer);
8651 sprintf(err_msg,
"BAD ATOM ID IN EXTRA BONDS FILE %s: %s",
8652 file->
data, buffer);
8660 int extraNumBonds = bond_params.
size();
8661 if ( extraNumBonds ) {
8662 iout <<
iINFO <<
"READ " << extraNumBonds <<
" EXTRA BONDS\n" <<
endi;
8664 #ifndef MEM_OPT_VERSION
8665 Bond *newbonds =
new Bond[numBonds+extraNumBonds];
8666 memcpy(newbonds, this->bonds, numBonds*
sizeof(
Bond));
8667 memcpy(newbonds+numBonds, bonds.
begin(), extraNumBonds*
sizeof(
Bond));
8668 delete [] this->bonds;
8669 this->bonds = newbonds;
8670 numBonds += extraNumBonds;
8684 int extraNumAngles = angle_params.
size();
8685 if ( extraNumAngles ) {
8686 iout <<
iINFO <<
"READ " << extraNumAngles <<
" EXTRA ANGLES\n" <<
endi;
8687 if ( anglesNormal ) {
8688 iout <<
iINFO <<
"EXTRA ANGLES ARE NORMAL HARMONIC\n" <<
endi;
8692 iout <<
iWARN <<
"EXTRA ANGLES ARE COSINE-BASED BY DEFAULT TO MATCH PREVIOUS VERSIONS\n";
8693 iout << iWARN <<
"FOR NORMAL HARMONIC EXTRA ANGLES SET extraBondsCosAngles off\n" <<
endi;
8695 #ifndef MEM_OPT_VERSION
8696 Angle *newangles =
new Angle[numAngles+extraNumAngles];
8697 memcpy(newangles, this->angles, numAngles*
sizeof(
Angle));
8698 memcpy(newangles+numAngles, angles.
begin(), extraNumAngles*
sizeof(
Angle));
8699 delete [] this->angles;
8700 this->angles = newangles;
8701 numAngles += extraNumAngles;
8715 int extraNumDihedrals = dihedral_params.
size();
8716 if ( extraNumDihedrals ) {
8717 iout <<
iINFO <<
"READ " << extraNumDihedrals <<
" EXTRA DIHEDRALS\n" <<
endi;
8718 #ifndef MEM_OPT_VERSION
8720 memcpy(newdihedrals, this->dihedrals, numDihedrals*
sizeof(
Dihedral));
8721 memcpy(newdihedrals+numDihedrals, dihedrals.
begin(), extraNumDihedrals*
sizeof(
Dihedral));
8722 delete [] this->dihedrals;
8723 this->dihedrals = newdihedrals;
8724 numDihedrals += extraNumDihedrals;
8738 int extraNumImpropers = improper_params.
size();
8739 if ( extraNumImpropers ) {
8740 iout <<
iINFO <<
"READ " << extraNumImpropers <<
" EXTRA IMPROPERS\n" <<
endi;
8741 #ifndef MEM_OPT_VERSION
8743 memcpy(newimpropers, this->impropers, numImpropers*
sizeof(
Improper));
8744 memcpy(newimpropers+numImpropers, impropers.
begin(), extraNumImpropers*
sizeof(
Improper));
8745 delete [] this->impropers;
8746 this->impropers = newimpropers;
8747 numImpropers += extraNumImpropers;
8781 PDB *initial_pdb,
char *cwd,
8782 const char *simmethod) {
8792 if (alchfile == NULL) {
8793 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, alchfile required.");
8795 strcpy(filename,
"coordinate pdb file (default)");
8798 if (alchfile->
next != NULL) {
8799 char *new_err_msg =
new char[24 + strlen(simmethod) + 26];
8800 sprintf(new_err_msg,
"Multiple definitions of %sFile in configuration file",simmethod);
8804 if ((cwd == NULL) || (alchfile->
data[0] ==
'/')) {
8805 strcpy(filename, alchfile->
data);
8808 strcpy(filename, cwd);
8809 strcat(filename, alchfile->
data);
8812 bPDB =
new PDB(filename);
8814 NAMD_die(
"Memory allocation failed in Molecule:build_fep_flags");
8818 char *new_err_msg =
new char[19 + strlen(simmethod) + 38];
8819 sprintf(new_err_msg,
"Number of atoms in %sFile PDB does not match coordinate PDB",simmethod);
8827 if (alchcol == NULL) {
8831 if (alchcol->
next != NULL) {
8832 char *new_err_msg =
new char[24 + strlen(simmethod) + 35];
8833 sprintf(new_err_msg,
"Multiple definitions of %s parameter column in config file",simmethod);
8837 if (strcasecmp(alchcol->
data,
"X") == 0) {
8840 else if (strcasecmp(alchcol->
data,
"Y") == 0) {
8843 else if (strcasecmp(alchcol->
data,
"Z") == 0) {
8846 else if (strcasecmp(alchcol->
data,
"O") == 0) {
8849 else if (strcasecmp(alchcol->
data,
"B") == 0) {
8853 NAMD_die(
"alchcol must have value of X, Y, Z, O or B");
8857 iout <<
iINFO <<
"To read " << simmethod <<
"data from file: " << filename
8859 iout <<
iINFO <<
"To read " << simmethod <<
"flag data from column: " << bcol
8863 fepAtomFlags =
new unsigned char[numAtoms];
8865 if (fepAtomFlags == NULL) {
8866 NAMD_die(
"Memory allocation failed in Molecule::build_fep_params()");
8869 double lesMassFactor = 1.0;
8871 lesMassFactor = 1.0 / simParams->
lesFactor;
8875 for (i = 0; i < numAtoms; i++) {
8879 bval = (bPDB->
atom(i))->xcoor();
8882 bval = (bPDB->
atom(i))->ycoor();
8885 bval = (bPDB->
atom(i))->zcoor();
8888 bval = (bPDB->
atom(i))->occupancy();
8891 bval = (bPDB->
atom(i))->temperaturefactor();
8896 if (simParams->
lesOn) {
8897 if ( bval == (
int) bval && bval > 0 ) {
8899 NAMD_die(
"LES flag must be less than or equal to lesFactor.");
8900 fepAtomFlags[i] = (int) bval;
8901 #ifdef MEM_OPT_VERSION
8902 Real newMass = atomMassPool[eachAtomMass[i]]*lesMassFactor;
8903 eachAtomMass[i] = insert_new_mass(newMass);
8905 atoms[i].mass *= lesMassFactor;
8910 fepAtomFlags[i] = 0;
8912 }
else if (simParams->
alchOn) {
8914 fepAtomFlags[i] = 3;
8916 }
else if (bval == 1.0) {
8919 }
else if (bval == -1.0) {
8920 fepAtomFlags[i] = 2;
8922 }
else if (bval == -2.0) {
8923 fepAtomFlags[i] = 4;
8926 fepAtomFlags[i] = 0;
8930 fepAtomFlags[i] = 1;
8933 fepAtomFlags[i] = 2;
8936 fepAtomFlags[i] = 0;
8939 fepAtomFlags[i] = (int) bval;
8941 #ifdef OPENATOM_VERSION
8943 if (simParams->openatomOn) {
8945 fepAtomFlags[i] = bval;
8948 fepAtomFlags[i] = 0;
8951 #endif //OPENATOM_VERSION
8955 if (alchfile != NULL) {
8979 if (ssfile == NULL) {
8980 if ( ! initial_pdb ) {
8981 NAMD_die(
"Initial PDB file unavailable, soluteScalingFile required.");
8984 strcpy(filename,
"coordinate PDB file (default)");
8987 if (ssfile->
next != NULL) {
8988 NAMD_die(
"Multiple definitions of soluteScalingFile in configuration file");
8991 if ((cwd == NULL) || (ssfile->
data[0] ==
'/')) {
8992 strcpy(filename, ssfile->
data);
8995 strcpy(filename, cwd);
8996 strcat(filename, ssfile->
data);
8999 bPDB =
new PDB(filename);
9001 NAMD_die(
"Memory allocation failed in Molecule::build_ss_flags");
9005 NAMD_die(
"Number of atoms in soluteScalingFile PDB does not match coordinate PDB");
9009 if (sscol == NULL) {
9013 if (sscol->
next != NULL) {
9014 NAMD_die(
"Multiple definitions of soluteScalingCol value in config file");
9017 if (strcasecmp(sscol->
data,
"X") == 0) {
9020 else if (strcasecmp(sscol->
data,
"Y") == 0) {
9023 else if (strcasecmp(sscol->
data,
"Z") == 0) {
9026 else if (strcasecmp(sscol->
data,
"O") == 0) {
9029 else if (strcasecmp(sscol->
data,
"B") == 0) {
9033 NAMD_die(
"soluteScalingCol must have value of X, Y, Z, O or B");
9037 iout <<
iINFO <<
"Reading solute scaling data from file: "
9038 << filename <<
"\n" <<
endi;
9039 iout <<
iINFO <<
"Reading solute scaling flags from column: "
9040 << bcol <<
"\n" <<
endi;
9042 ssAtomFlags =
new unsigned char[numAtoms];
9043 ss_index =
new int[numAtoms];
9045 if (ssAtomFlags == NULL || ss_index == NULL) {
9046 NAMD_die(
"Memory allocation failed in Molecule::build_ss_params()");
9050 for (i = 0; i < numAtoms; i++) {
9053 bval = (bPDB->
atom(i))->xcoor();
9056 bval = (bPDB->
atom(i))->ycoor();
9059 bval = (bPDB->
atom(i))->zcoor();
9062 bval = (bPDB->
atom(i))->occupancy();
9065 bval = (bPDB->
atom(i))->temperaturefactor();
9071 ss_index[num_ss] = i;
9080 if (ssfile != NULL) {
9089 int *numAtomsByLjType =
new int[LJtypecount];
9093 ss_vdw_type =
new int[LJtypecount];
9096 for (i = 0; i < LJtypecount; i++) {
9097 numAtomsByLjType[i] = 0;
9102 for (i = 0; i < num_ss; i++) {
9103 numAtomsByLjType[atoms[ss_index[i]].vdw_type]++;
9107 ss_num_vdw_params = 0;
9108 for (i = 0; i < LJtypecount; i++) {
9110 if (numAtomsByLjType[i] != 0) {
9113 ss_vdw_type[ss_num_vdw_params] = i;
9116 ss_num_vdw_params++;
9120 for (i = 0; i < num_ss; i++) {
9122 for (j = 0; j < ss_num_vdw_params; j++) {
9124 if (atoms[ss_index[i]].vdw_type == ss_vdw_type[j]) {
9127 atoms[ss_index[i]].vdw_type = LJtypecount + j;
9132 delete [] numAtomsByLjType;
9146 #ifndef MEM_OPT_VERSION
9150 FILE *alch_unpert_bond_file;
9153 if ((alch_unpert_bond_file =
Fopen(alch_fname,
"r")) == NULL) {
9154 sprintf(err_msg,
"UNABLE TO OPEN ALCH UNPERTBURBED BOND FILE %s", alch_fname);
9164 NAMD_die(
"DID NOT FIND NBOND AFTER ATOM LIST IN ALCH UNPERT PSF");
9168 sscanf(buffer,
"%d", &num_alch_unpert_Bonds);
9170 read_alch_unpert_bonds(alch_unpert_bond_file);
9179 NAMD_die(
"DID NOT FIND NTHETA AFTER BOND LIST IN ALCH UNPERT PSF");
9183 sscanf(buffer,
"%d", &num_alch_unpert_Angles);
9185 read_alch_unpert_angles(alch_unpert_bond_file);
9196 NAMD_die(
"DID NOT FIND NPHI AFTER ANGLE LIST IN ALCH UNPERT PSF");
9200 sscanf(buffer,
"%d", &num_alch_unpert_Dihedrals);
9202 read_alch_unpert_dihedrals(alch_unpert_bond_file);
9204 Fclose(alch_unpert_bond_file);
9210 suspiciousAlchBonds = 0;
9211 for (
int i = 0; i < numBonds; i++) {
9212 int part1 = fepAtomFlags[bonds[i].atom1];
9213 int part2 = fepAtomFlags[bonds[i].atom2];
9214 if ((part1 == 1 || part2 == 1 ) &&
9215 (part1 == 2 || part2 == 2 )) {
9217 suspiciousAlchBonds++;
9222 Angle *nonalchAngles;
9223 nonalchAngles =
new Angle[numAngles];
9224 int nonalchAngleCount = 0;
9225 alchDroppedAngles = 0;
9226 for (
int i = 0; i < numAngles; i++) {
9227 int part1 = fepAtomFlags[angles[i].atom1];
9228 int part2 = fepAtomFlags[angles[i].atom2];
9229 int part3 = fepAtomFlags[angles[i].atom3];
9230 if ((part1 == 1 || part2 == 1 || part3 == 1) &&
9231 (part1 == 2 || part2 == 2 || part3 == 2)) {
9233 alchDroppedAngles++;
9236 if ( angles[i].angle_type == -1 ) {
9239 "MISSING PARAMETERS FOR ANGLE %i %i %i PARTITIONS %i %i %i\n",
9240 angles[i].atom1+1, angles[i].atom2+1, angles[i].atom3+1,
9241 part1, part2, part3);
9244 nonalchAngles[nonalchAngleCount++] = angles[i];
9247 numAngles = nonalchAngleCount;
9249 angles =
new Angle[numAngles];
9250 for (
int i = 0; i < nonalchAngleCount; i++) {
9251 angles[i]=nonalchAngles[i];
9253 delete [] nonalchAngles;
9258 nonalchDihedrals =
new Dihedral[numDihedrals];
9259 int nonalchDihedralCount = 0;
9260 alchDroppedDihedrals = 0;
9261 for (
int i = 0; i < numDihedrals; i++) {
9262 int part1 = fepAtomFlags[dihedrals[i].atom1];
9263 int part2 = fepAtomFlags[dihedrals[i].atom2];
9264 int part3 = fepAtomFlags[dihedrals[i].atom3];
9265 int part4 = fepAtomFlags[dihedrals[i].atom4];
9266 if ((part1 == 1 || part2 == 1 || part3 == 1 || part4 == 1) &&
9267 (part1 == 2 || part2 == 2 || part3 == 2 || part4 == 2)) {
9269 alchDroppedDihedrals++;
9272 if ( dihedrals[i].dihedral_type == -1 ) {
9275 "MISSING PARAMETERS FOR DIHEDRAL %i %i %i %i PARTITIONS %i %i %i %i\n",
9276 dihedrals[i].atom1+1, dihedrals[i].atom2+1,
9277 dihedrals[i].atom3+1, dihedrals[i].atom4+1,
9278 part1, part2, part3, part4);
9281 nonalchDihedrals[nonalchDihedralCount++] = dihedrals[i];
9284 numDihedrals = nonalchDihedralCount;
9285 delete [] dihedrals;
9286 dihedrals =
new Dihedral[numDihedrals];
9287 for (
int i = 0; i < numDihedrals; i++) {
9288 dihedrals[i]=nonalchDihedrals[i];
9290 delete [] nonalchDihedrals;
9294 nonalchImpropers =
new Improper[numImpropers];
9295 int nonalchImproperCount = 0;
9296 alchDroppedImpropers = 0;
9297 for (
int i = 0; i < numImpropers; i++) {
9298 int part1 = fepAtomFlags[impropers[i].atom1];
9299 int part2 = fepAtomFlags[impropers[i].atom2];
9300 int part3 = fepAtomFlags[impropers[i].atom3];
9301 int part4 = fepAtomFlags[impropers[i].atom4];
9302 if ((part1 == 1 || part2 == 1 || part3 == 1 || part4 == 1) &&
9303 (part1 == 2 || part2 == 2 || part3 == 2 || part4 == 2)) {
9305 alchDroppedImpropers++;
9308 nonalchImpropers[nonalchImproperCount++] = impropers[i];
9311 numImpropers = nonalchImproperCount;
9312 delete [] impropers;
9313 impropers =
new Improper[numImpropers];
9314 for (
int i = 0; i < numImpropers; i++) {
9315 impropers[i]=nonalchImpropers[i];
9317 delete [] nonalchImpropers;
9338 if (fixedfile == NULL) {
9339 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, excludeFromPressureFile required.");
9342 if (fixedfile->
next != NULL) {
9343 NAMD_die(
"Multiple definitions of excluded pressure atoms PDB file in configuration file");
9346 if ( (cwd == NULL) || (fixedfile->
data[0] ==
'/') ) {
9347 strcpy(filename, fixedfile->
data);
9349 strcpy(filename, cwd);
9350 strcat(filename, fixedfile->
data);
9352 bPDB =
new PDB(filename);
9353 if ( bPDB == NULL ) {
9354 NAMD_die(
"Memory allocation failed in Molecule::build_exPressure_atoms");
9358 NAMD_die(
"Number of atoms in excludedPressure atoms PDB doesn't match coordinate PDB");
9367 if (fixedcol == NULL) {
9370 if (fixedcol->
next != NULL) {
9371 NAMD_die(
"Multiple definitions of excludedPressure atoms column in config file");
9374 if (strcasecmp(fixedcol->
data,
"X") == 0) {
9376 }
else if (strcasecmp(fixedcol->
data,
"Y") == 0) {
9378 }
else if (strcasecmp(fixedcol->
data,
"Z") == 0) {
9380 }
else if (strcasecmp(fixedcol->
data,
"O") == 0) {
9382 }
else if (strcasecmp(fixedcol->
data,
"B") == 0) {
9385 NAMD_die(
"excludedPressureFileCol must have value of X, Y, Z, O, or B");
9390 exPressureAtomFlags =
new int32[numAtoms];
9392 if (exPressureAtomFlags == NULL) {
9393 NAMD_die(
"memory allocation failed in Molecule::build_fixed_atoms()");
9396 numExPressureAtoms = 0;
9399 for (i=0; i<numAtoms; i++) {
9402 case 1: bval = (bPDB->
atom(i))->xcoor();
break;
9403 case 2: bval = (bPDB->
atom(i))->ycoor();
break;
9404 case 3: bval = (bPDB->
atom(i))->zcoor();
break;
9405 case 4: bval = (bPDB->
atom(i))->occupancy();
break;
9406 case 5: bval = (bPDB->
atom(i))->temperaturefactor();
break;
9411 exPressureAtomFlags[i] = 1;
9412 numExPressureAtoms++;
9414 exPressureAtomFlags[i] = 0;
9417 if (fixedfile != NULL)
9420 iout <<
iINFO <<
"Got " << numExPressureAtoms <<
" excluded pressure atoms."
9426 return ((atoms[anum].status & LonepairAtom) != 0);
9430 return ((atoms[anum].status &
DrudeAtom) != 0);
9440 return ((atoms[anum].status &
OxygenAtom) != 0);
9445 return (hydrogenGroup[atoms[anum].hydrogenList].isGP);
9450 return (hydrogenGroup[atoms[anum].hydrogenList].waterVal == 2);
9455 return (hydrogenGroup[atoms[anum].hydrogenList].atomsInGroup);
9462 return atoms[anum].partner;
9466 if ( n != numAtoms )
9467 NAMD_die(
"Incorrect number of atoms in Molecule::reloadCharges().");
9469 #ifdef MEM_OPT_VERSION
9470 delete [] atomChargePool;
9471 vector<Real> tmpCharges;
9472 for(
int i=0; i<numAtoms; i++){
9476 for(
int j=0; j<tmpCharges.size();j++){
9477 if(tmpCharges[j] == charge[i]){
9483 tmpCharges.push_back(charge[i]);
9484 foundIdx = tmpCharges.size()-1;
9486 eachAtomCharge[i] = (
Index)foundIdx;
9488 chargePoolSize = tmpCharges.size();
9489 atomChargePool =
new Real[chargePoolSize];
9490 for(
int i=0; i<chargePoolSize; i++)
9491 atomChargePool[i] = tmpCharges[i];
9493 for(
int i=0; i<n; ++i ) atoms[i].charge = charge[i];
9497 #ifndef MEM_OPT_VERSION
9504 void Molecule::build_atom_status(
void) {
9507 int numDrudeWaters = 0;
9510 numLonepairs = numZeroMassAtoms;
9512 iout <<
iWARN <<
"CORRECTION OF ZERO MASS ATOMS TURNED OFF "
9513 "BECAUSE LONE PAIRS ARE USED\n" <<
endi;
9517 if (is_lonepairs_psf && numLonepairs != numLphosts) {
9518 NAMD_die(
"must have same number of LP hosts as lone pairs");
9520 }
else if (numZeroMassAtoms) {
9521 for (i=0; i < numAtoms; i++) {
9522 if ( atoms[i].mass < 0.001 ) atoms[i].mass = 0.001;
9525 iout <<
iWARN <<
"FOUND " << numZeroMassAtoms <<
9526 " ATOMS WITH ZERO OR NEGATIVE MASSES! CHANGED TO 0.001\n" <<
endi;
9531 hydrogenGroup.resize(numAtoms);
9533 for (i=0; i < numAtoms; i++) {
9534 atoms[i].partner = (-1);
9544 int hhbondcount = 0;
9545 for (i=0; i < numRealBonds; i++) {
9546 a1 = bonds[i].atom1;
9547 a2 = bonds[i].atom2;
9548 if (is_hydrogen(a1) && is_hydrogen(a2)) {
9551 atoms[a1].partner = a2;
9552 atoms[a2].partner = a1;
9560 if ( hhbondcount && ! CkMyPe() ) {
9561 iout <<
iWARN <<
"Found " << hhbondcount <<
" H-H bonds.\n" <<
endi;
9566 for (i=0; i < numRealBonds; i++) {
9567 a1 = bonds[i].atom1;
9568 a2 = bonds[i].atom2;
9569 if (is_hydrogen(a1)) {
9570 if (is_hydrogen(a2))
continue;
9571 atoms[a1].partner = a2;
9577 if (is_oxygen(a2)) hg[a2].waterVal++;
9579 if (is_hydrogen(a2)) {
9580 atoms[a2].partner = a1;
9586 if (is_oxygen(a1)) hg[a1].waterVal++;
9592 atoms[a1].partner = a2;
9599 atoms[a2].partner = a1;
9607 else if ( is_lonepairs_psf || is_drude_psf ) {
9608 if (is_lp(a1) || is_drude(a1)) {
9609 if (is_hydrogen(a2) || is_lp(a2) || is_drude(a2)) {
9611 sprintf(msg,
"%s particle %d is bonded to non-parent atom %d",
9612 (is_lp(a1) ?
"Lone pair" :
"Drude"), a1+1, a2+1);
9615 atoms[a1].partner = a2;
9621 else if (is_lp(a2) || is_drude(a2)) {
9622 if (is_hydrogen(a1) || is_lp(a1) || is_drude(a1)) {
9624 sprintf(msg,
"%s particle %d is bonded to non-parent atom %d",
9625 (is_lp(a2) ?
"Lone pair" :
"Drude"), a2+1, a1+1);
9628 atoms[a2].partner = a1;
9640 for(i=0; i<numAtoms; i++) {
9641 if ( ! hg[hg[i].GPID].isGP ) {
9643 sprintf(msg,
"child atom %d bonded only to child H atoms",i+1);
9646 if ( hg[i].isGP && is_hydrogen(i) ) {
9647 if ( hg[i].GPID == i )
continue;
9649 if ( is_hydrogen(hg[i].GPID) && hg[hg[i].GPID].GPID != i ) {
9651 sprintf(msg,
"H atom %d bonded only to child H atoms",i+1);
9657 if ( hg[i].atomsInGroup != 2 ) {
9659 sprintf(msg,
"H atom %d bonded to multiple H atoms",i+1);
9664 if ( hGPcount && ! CkMyPe() ) {
9665 iout <<
iWARN <<
"Found " << hGPcount <<
" H-H molecules.\n" <<
endi;
9669 for (i=0; i<numAtoms; ++i) {
9670 if ( hg[i].isGP ) hg[i].
GPID = i;
9678 for (i=0; i<numLphosts; ++i) {
9679 int a1 = lphosts[i].atom1;
9680 int a2 = lphosts[i].atom2;
9681 int a3 = lphosts[i].atom3;
9682 int a4 = lphosts[i].atom4;
9683 int m1 = hg[a1].
MPID;
9684 while ( hg[m1].MPID != m1 ) m1 = hg[m1].
MPID;
9685 int m2 = hg[a2].
MPID;
9686 while ( hg[m2].MPID != m2 ) m2 = hg[m2].
MPID;
9687 int m3 = hg[a3].
MPID;
9688 while ( hg[m3].MPID != m3 ) m3 = hg[m3].
MPID;
9689 int m4 = hg[a4].
MPID;
9690 while ( hg[m4].MPID != m4 ) m4 = hg[m4].
MPID;
9692 if ( m2 < mp ) mp = m2;
9693 if ( m3 < mp ) mp = m3;
9694 if ( m4 < mp ) mp = m4;
9702 for (i=0; i<numAtoms; ++i) {
9703 int mp = hg[i].
MPID;
9704 if ( hg[mp].MPID != mp ) {
9711 for (i=0; i<numAtoms; ++i) {
9715 for (i=0; i<numAtoms; ++i) {
9721 for (i=0; i<numAtoms; i++) {
9732 numHydrogenGroups = 0;
9733 maxHydrogenGroupSize = 0;
9734 numMigrationGroups = 0;
9735 maxMigrationGroupSize = 0;
9736 for(i=0; i<numAtoms; i++)
9739 ++numMigrationGroups;
9741 if ( mgs > maxMigrationGroupSize ) maxMigrationGroupSize = mgs;
9744 ++numHydrogenGroups;
9746 if ( hgs > maxHydrogenGroupSize ) maxHydrogenGroupSize = hgs;
9750 hydrogenGroup.sort();
9755 for(i=0; i<numAtoms; ++i, --hgs) {
9762 sprintf(buff,
"Atom %d has bad hydrogen group size. "
9763 "Check for duplicate bonds.", parentid+1);
9769 sprintf(buff,
"Atom %d has bad hydrogen group size. "
9770 "Check for duplicate bonds.", parentid+1);
9778 for(i=0; i<numAtoms; ++i, --mgs) {
9785 sprintf(buff,
"Atom %d has bad migration group size.", parentid+1);
9791 sprintf(buff,
"Atom %d has bad migration group size.", parentid+1);
9799 for(i=0; i<numAtoms; i++) {
9800 atoms[hydrogenGroup[i].atomID].hydrogenList = i;
9806 for (i = 0; i < numAtoms; i++) {
9807 if (is_water(hg[i].atomID) && hg[i].isGP) {
9809 || ! is_drude(hg[i+1].atomID)
9810 || ! is_lp(hg[i+2].atomID)
9811 || ! is_hydrogen(hg[i+3].atomID)
9812 || ! is_hydrogen(hg[i+4].atomID) ) {
9814 sprintf(msg,
"Drude water molecule from HydrogenGroup i=%d "
9815 "starting at atom %d is not sorted\n", i, hg[i].atomID+1);
9822 else if (is_drude(hg[i].atomID)) {
9823 if (i < 1 || hg[i-1].atomID != hg[i].GPID) {
9825 sprintf(msg,
"Drude particle from HydrogenGroup i=%d must "
9826 "immediately follow its parent atom %d\n", i, hg[i].GPID+1);
9831 else if (is_lp(hg[i].atomID)) {
9833 sprintf(msg,
"Drude lonepair from HydrogenGroup i=%d "
9834 "at particle %d is NOT from water - unsupported\n",
9845 for(i=0; i<numAtoms; i++)
9846 iout << i <<
" atomID=" << hydrogenGroup[i].atomID
9847 <<
" isGP=" << hydrogenGroup[i].isGP
9848 <<
" parent=" << hydrogenGroup[i].GPID
9849 <<
" #" << hydrogenGroup[i].atomsInGroup
9850 <<
" waterVal=" << hydrogenGroup[i].waterVal
9851 <<
" partner=" << atoms[i].partner
9852 <<
" hydrogenList=" << atoms[i].hydrogenList
9863 delete [] rigidBondLengths;
9864 rigidBondLengths =
new Real[numAtoms];
9865 if ( ! rigidBondLengths ) {
9866 NAMD_die(
"Memory allocation failed in Molecule::build_atom_status()\n");
9868 for (i=0; i<numAtoms; ++i) rigidBondLengths[i] = 0;
9873 for (i=0; i < numRealBonds; i++) {
9874 a1 = bonds[i].atom1;
9875 a2 = bonds[i].atom2;
9878 if (is_hydrogen(a2)) {
int tmp = a1; a1 = a2; a2 = tmp; }
9879 if (is_hydrogen(a1)) {
9880 if ( is_hydrogen(a2) ) {
9881 if ( ! is_water(a2) ) {
9882 rigidBondLengths[a1] = ( mode ==
RIGID_ALL ? x0 : 0. );
9883 rigidBondLengths[a2] = ( mode ==
RIGID_ALL ? x0 : 0. );
9885 }
else if ( is_water(a2) || mode ==
RIGID_ALL ) {
9886 rigidBondLengths[a1] = x0;
9887 if (is_water(a2)) r_oh = rigidBondLengths[a1];
9889 rigidBondLengths[a1] = 0.;
9894 if (is_lp(a2)) {
int tmp = a1; a1 = a2; a2 = tmp; }
9896 if (! is_water(a2) ) {
9900 sprintf(err_msg,
"ILLEGAL LONE PAIR AT INDEX %i\n"
9901 "LONE PAIRS ARE CURRENTLY ALLOWED ONLY ON WATER MOLECULES\n",
9905 rigidBondLengths[a1] = x0;
9914 int tmp = a1; a1 = a2; a2 = tmp;
9921 else if ( ! simParams->
drudeOn) {
9924 sprintf(msg,
"ILLEGAL LONE PAIR AT INDEX %d\n"
9925 "LONE PAIRS ARE CURRENTLY ALLOWED ONLY ON WATER MOLECULES\n",
9935 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
9936 for( ; h_i != h_e; ++h_i ) {
9937 if ( h_i->isGP ) rigidBondLengths[h_i->atomID] = 0.;
9941 for (i=0; i < numAngles; i++) {
9942 a2 = angles[i].atom2;
9943 if ( ! is_water(a2) )
continue;
9944 if ( ! is_oxygen(a2) )
continue;
9945 a1 = angles[i].atom1;
9946 if ( ! is_hydrogen(a1) )
continue;
9947 a3 = angles[i].atom3;
9948 if ( ! is_hydrogen(a3) )
continue;
9949 if (is_lp(a2) || is_lp(a1) || is_lp(a3) ||
9950 is_drude(a2) || is_drude(a1) || is_drude(a3))
continue;
9951 if ( rigidBondLengths[a1] != rigidBondLengths[a3] ) {
9952 if (rigidBondLengths[a1] >0.3 && rigidBondLengths[a3] >0.3) {
9953 printf(
"length1: %f length2: %f\n", rigidBondLengths[a1], rigidBondLengths[a3]);
9955 NAMD_die(
"Asymmetric water molecule found??? This can't be right.\n");
9960 rigidBondLengths[a2] = 2. * rigidBondLengths[a1] * sin(0.5*t0);
9961 r_hh = rigidBondLengths[a2];
9965 int numBondWaters = 0;
9966 int numFailedWaters = 0;
9968 for (i=0; i < numRealBonds; i++) {
9969 a1 = bonds[i].atom1;
9970 a2 = bonds[i].atom2;
9971 if ( ! is_hydrogen(a1) )
continue;
9972 if ( ! is_hydrogen(a2) )
continue;
9973 int ma1 = get_mother_atom(a1);
9974 int ma2 = get_mother_atom(a2);
9975 if ( ma1 != ma2 )
continue;
9976 if ( ! is_water(ma1) )
continue;
9977 if ( rigidBondLengths[ma1] != 0. )
continue;
9980 rigidBondLengths[ma1] = x0;
9987 if (r_oh < 0.0 || r_hh < 0.0) {
9989 NAMD_die(
"Failed to find water bond lengths\n");
9991 r_ohc = sqrt(r_oh * r_oh - 0.25 * r_hh * r_hh);
9995 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
9996 for( ; h_i != h_e; ++h_i ) {
9997 if ( h_i->isGP && is_water(h_i->atomID) &&
9998 rigidBondLengths[h_i->atomID] == 0. ) {
9999 if ( h_i + 1 == h_e || h_i + 2 == h_e ||
10000 h_i[1].isGP || h_i[2].isGP || h_i->atomsInGroup != 3 ) {
10001 NAMD_die(
"Abnormal water detected.");
10003 if ( CkNumNodes() > 1 ) {
10004 NAMD_die(
"Unable to determine H-H distance for rigid water because structure has neither H-O-H angle nor H-H bond.");
10007 btmp.
atom1 = h_i[1].atomID;
10008 btmp.
atom2 = h_i[2].atomID;
10010 get_atomtype(btmp.
atom1),
10011 get_atomtype(btmp.
atom2),
10017 rigidBondLengths[h_i->atomID] = x0;
10024 if ( numBondWaters + numFailedWaters ) {
10025 iout <<
iWARN <<
"Missing angles for " <<
10026 ( numBondWaters + numFailedWaters ) <<
" waters.\n" << endi;
10028 if ( numBondWaters ) {
10029 iout <<
iWARN <<
"Obtained H-H distance from bond parameters for " <<
10030 numBondWaters <<
" waters.\n" <<
endi;
10031 iout <<
iWARN <<
"This would not be possible in a multi-process run.\n" <<
endi;
10033 if ( numFailedWaters ) {
10034 iout <<
iERROR <<
"Failed to obtain H-H distance from angles or bonds for " <<
10035 numFailedWaters <<
" waters.\n" <<
endi;
10043 for (i=0; i<numAtoms; ++i) rigidBondLengths[i] *= -1;
10045 for (i=0; i<numAtoms; ++i) {
10046 if ( ! is_water(i) ) rigidBondLengths[i] *= -1;
10052 for (i=0; i<numAtoms; ++i) {
10053 if ( rigidBondLengths[i] > 0. ) ++numRigidBonds;
10074 BigReal LJAvgA, LJAvgB, LJAvgA1, LJAvgB1, LJAvgA2, LJAvgB2;
10075 int numLJsites, numLJsites1, numLJsites2;
10085 Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
10086 Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
10087 Real *ATable =
new Real[LJtypecount*LJtypecount];
10088 Real *BTable =
new Real[LJtypecount*LJtypecount];
10091 for (
int i = 0; i < LJtypecount; i++) {
10092 for (
int j = 0; j < LJtypecount; j++) {
10094 ATable[i*LJtypecount + j] =
A;
10095 BTable[i*LJtypecount + j] =
B;
10098 params->
get_vdw_params(&sigma_i,&epsilon_i,&sigma_i14,&epsilon_i14,i);
10099 params->
get_vdw_params(&sigma_j,&epsilon_j,&sigma_j14,&epsilon_j14,j);
10101 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10102 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10103 sigma_ij *= sigma_ij*sigma_ij;
10104 sigma_ij *= sigma_ij;
10106 ATable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10107 BTable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij;
10112 int *numAtomsByLjType =
new int[LJtypecount];
10113 for (
int i=0; i < LJtypecount; i++) {numAtomsByLjType[i]=0;}
10114 for (
int i=0; i < numAtoms; i++) {numAtomsByLjType[atoms[i].vdw_type]++;}
10121 for (
int i=0; i < LJtypecount; i++) {
10122 for (
int j=0; j < LJtypecount; j++) {
10123 A = ATable[i*LJtypecount + j];
10124 B = BTable[i*LJtypecount + j];
10125 if (!A && !B)
continue;
10126 npairs = (numAtomsByLjType[i] - int(i==j))*
BigReal(numAtomsByLjType[j]);
10127 sumOfAs += npairs*
A;
10128 sumOfBs += npairs*
B;
10130 if (i==j) numLJsites += numAtomsByLjType[i];
10133 delete [] numAtomsByLjType;
10137 LJAvgA = sumOfAs / count;
10138 LJAvgB = sumOfBs / count;
10159 numLJsites1 = numLJsites2 = numLJsites;
10160 int alch_counter = 0;
10161 for (
int i=0; i < numAtoms; ++i) {
10163 if (get_fep_type(i) == 2 || get_fep_type(i) == 4) alchFlagi = -1;
10164 if (get_fep_type(i) == 1 || get_fep_type(i) == 3) alchFlagi = 1;
10166 &A, &B, &A14, &B14)) {
10170 &epsilon_i14, atoms[i].vdw_type);
10172 useGeom ? sqrt(sigma_i*sigma_i) : 0.5*(sigma_i+sigma_i);
10173 BigReal epsilon_ii = sqrt(epsilon_i*epsilon_i);
10175 sigma_ii *= sigma_ii*sigma_ii;
10176 sigma_ii *= sigma_ii;
10177 A = 4.0*sigma_ii*epsilon_ii*sigma_ii;
10178 B = 4.0*sigma_ii*epsilon_ii;
10181 if (alchFlagi == 1) numLJsites2--;
10182 else if (alchFlagi == -1) numLJsites1--;
10184 for (
int j=i+1; j < numAtoms; ++j) {
10186 if (get_fep_type(j) == 2 || get_fep_type(j) == 4) alchFlagj = -1;
10187 if (get_fep_type(j) == 1 || get_fep_type(j) == 3) alchFlagj = 1;
10188 int alchFlagSum = alchFlagi + alchFlagj;
10191 if (alchFlagi == 0 && alchFlagj == 0)
continue;
10194 &A, &B, &A14, &B14)) {
10198 &epsilon_i14, atoms[i].vdw_type);
10200 &epsilon_j14, atoms[j].vdw_type);
10202 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10203 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10205 sigma_ij *= sigma_ij*sigma_ij;
10206 sigma_ij *= sigma_ij;
10207 A = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10208 B = 4.0*sigma_ij*epsilon_ij;
10210 if (!A && !B)
continue;
10215 if ( alchFlagSum > 0 ){
10220 else if ( alchFlagSum < 0 ){
10236 if ( alchFlagi == 1 || alchFlagi == -1 ) alch_counter++;
10237 if ( alch_counter == (numFepInitial + numFepFinal) )
break;
10239 LJAvgA = sumOfAs / count;
10240 LJAvgB = sumOfBs / count;
10242 LJAvgA1 = sumOfAs1 / count1;
10243 LJAvgB1 = sumOfBs1 / count1;
10249 LJAvgA2 = sumOfAs2 / count2;
10250 LJAvgB2 = sumOfBs2 / count2;
10255 if ( ! CkMyPe() ) {
10256 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO "
10257 <<
"ENERGY AND PRESSURE\n" <<
endi;
10258 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A0 AND B0 COEFFICIENTS "
10259 << LJAvgA2 <<
" AND " << LJAvgB2 <<
"\n" <<
endi;
10260 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A1 AND B1 COEFFICIENTS "
10261 << LJAvgA1 <<
" AND " << LJAvgB1 <<
"\n" <<
endi;
10263 numLJsites = (numLJsites1 + numLJsites2 - numLJsites);
10264 LJAvgA1 *=
BigReal(numLJsites1)*numLJsites1;
10265 LJAvgB1 *=
BigReal(numLJsites1)*numLJsites1;
10266 LJAvgA2 *=
BigReal(numLJsites2)*numLJsites2;
10267 LJAvgB2 *=
BigReal(numLJsites2)*numLJsites2;
10270 LJAvgA1 = LJAvgB1 = LJAvgA2 = LJAvgB2 = 0;
10272 if ( ! CkMyPe() ) {
10273 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO "
10274 <<
"ENERGY AND PRESSURE\n" <<
endi;
10275 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A AND B COEFFICIENTS "
10276 << LJAvgA <<
" AND " << LJAvgB <<
"\n" <<
endi;
10279 LJAvgA *=
BigReal(numLJsites)*numLJsites;
10280 LJAvgB *=
BigReal(numLJsites)*numLJsites;
10289 BigReal rswitch2 = rswitch*rswitch;
10290 BigReal rswitch3 = rswitch*rswitch2;
10291 BigReal rswitch4 = rswitch2*rswitch2;
10292 BigReal rswitch5 = rswitch2*rswitch3;
10293 BigReal rswitch6 = rswitch3*rswitch3;
10319 BigReal int_U_gofr_A, int_rF_gofr_A, int_U_gofr_B, int_rF_gofr_B;
10338 BigReal rsum3 = (rcut + rswitch)*(rcut + rswitch)*(rcut + rswitch);
10339 int_U_gofr_A = int_rF_gofr_A = (16*
PI*(3*rcut4 + 9*rcut3*rswitch
10340 + 11*rcut2*rswitch2 + 9*rcut*rswitch3
10342 / (315*rcut5*rswitch5*rsum3));
10343 int_U_gofr_B = int_rF_gofr_B = -16*
PI / (3*rsum3);
10370 BigReal lnr = log(rcut/rswitch);
10394 int_U_gofr_A = int_rF_gofr_A =\
10395 16*
PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
10396 int_U_gofr_B = int_rF_gofr_B = -4*
PI*lnr / (rcut3 - rswitch3);
10433 int_rF_gofr_A = 8*
PI / (9*rcut9);
10434 int_rF_gofr_B = -4*
PI / (3*rcut3);
10435 int_U_gofr_A = 2*
PI / (9*rcut9);
10436 int_U_gofr_B = -2*
PI / (3*rcut3);
10439 tail_corr_virial = int_rF_gofr_A*LJAvgA + int_rF_gofr_B*LJAvgB;
10440 tail_corr_ener = int_U_gofr_A*LJAvgA + int_U_gofr_B*LJAvgB;
10442 tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1 -
10444 tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1 -
10446 tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2 -
10448 tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2 -
10460 const BigReal corr = (doti ? 0.0 : tail_corr_ener);
10461 return scl*(corr + vdw_lambda_1*tail_corr_dUdl_1 +
10462 vdw_lambda_2*tail_corr_dUdl_2);
10465 return scl*tail_corr_ener;
10475 return scl*(tail_corr_virial + vdw_lambda_1*tail_corr_virial_1 +
10476 vdw_lambda_2*tail_corr_virial_2);
10479 return scl*tail_corr_virial;
10484 #ifdef MEM_OPT_VERSION
10487 int Molecule::checkExclByIdx(
int idx1,
int atom1,
int atom2)
const {
10489 int amin = exclChkSigPool[idx1].min;
10490 int amax = exclChkSigPool[idx1].max;
10491 int dist21 = atom2 - atom1;
10492 if ( dist21 < amin || dist21 > amax )
return 0;
10493 else return exclChkSigPool[idx1].flags[dist21-amin];
10499 int amin = all_exclusions[atom1].min;
10500 int amax = all_exclusions[atom1].max;
10501 if ( atom2 < amin || atom2 > amax )
return 0;
10502 else return all_exclusions[atom1].flags[atom2-amin];
10518 initialize(simParams,param);
10522 is_lonepairs_psf = 0;
10530 read_parm(amber_data);
10532 #ifndef MEM_OPT_VERSION
10535 assignLCPOTypes( 1 );
10543 initialize(simParams, param);
10547 is_lonepairs_psf = 0;
10555 read_parm(amber_data);
10557 #ifndef MEM_OPT_VERSION
10560 assignLCPOTypes( 1 );
10565 #ifdef MEM_OPT_VERSION
10566 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
10568 int i,j,ntheth,nphih,current_index,a1,a2,
10569 max,min,index,found;
10572 NAMD_die(
"No data read from parm file yet!");
10575 numAtoms = amber_data->
Natom;
10576 atoms =
new Atom[numAtoms];
10588 if (atoms == NULL || atomNames == NULL )
10589 NAMD_die(
"memory allocation failed when reading atom information");
10592 for (i = 0; i < numAtoms; ++i) {
10593 const int resname_length = amber_data->
ResNames[amber_data->
AtomRes[i]].size();
10594 const int atomname_length = amber_data->
AtomNames[i].size();
10595 const int atomtype_length = amber_data->
AtomSym[i].size();
10596 atomNames[i].resname = nameArena->getNewArray(resname_length+1);
10597 atomNames[i].atomname = nameArena->getNewArray(atomname_length+1);
10598 atomNames[i].atomtype = nameArena->getNewArray(atomtype_length+1);
10599 if (atomNames[i].resname == NULL || atomNames[i].atomname == NULL || atomNames[i].atomtype == NULL)
10600 NAMD_die(
"memory allocation failed when reading atom information");
10602 strncpy(atomNames[i].resname, amber_data->
ResNames[amber_data->
AtomRes[i]].c_str(), resname_length+1);
10603 strncpy(atomNames[i].atomname, amber_data->
AtomNames[i].c_str(), atomname_length+1);
10604 strncpy(atomNames[i].atomtype, amber_data->
AtomSym[i].c_str(), atomtype_length+1);
10606 strtok(atomNames[i].resname,
" ");
10607 strtok(atomNames[i].atomname,
" ");
10608 strtok(atomNames[i].atomtype,
" ");
10609 atoms[i].mass = amber_data->
Masses[i];
10656 atoms[i].charge = amber_data->
Charges[i] / 18.2223;
10657 atoms[i].vdw_type = amber_data->
Iac[i] - 1;
10659 if ( tmpResLookup ) tmpResLookup =
10662 if(atomSegResids) {
10664 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
10671 }
else if (atoms[i].mass <= 0.05) {
10673 ++numZeroMassAtoms;
10674 }
else if (atoms[i].mass < 1.0) {
10680 }
else if (atoms[i].mass <=3.5) {
10682 }
else if ((atomNames[i].atomname[0] ==
'O') &&
10683 (atoms[i].mass >= 14.0) &&
10684 (atoms[i].mass <= 18.0)) {
10697 if (amber_data->
Nbonh + amber_data->
Nbona > 0) {
10699 if (bonds == NULL || amber_data->
Nbona < 0)
10700 NAMD_die(
"memory allocation failed when reading bond information");
10702 for (i=0; i<amber_data->
Nbonh; ++i)
10703 { bonds[numBonds].atom1 = amber_data->
BondHAt1[i] / 3;
10704 bonds[numBonds].atom2 = amber_data->
BondHAt2[i] / 3;
10705 bonds[numBonds].bond_type = amber_data->
BondHNum[i] - 1;
10706 if (bonds[numBonds].atom1>=numAtoms || bonds[numBonds].atom2>=numAtoms ||
10707 bonds[numBonds].bond_type>=amber_data->
Numbnd)
10708 {
char err_msg[128];
10709 sprintf(err_msg,
"BOND (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
10717 for (i=amber_data->
Nbonh; i<amber_data->Nbonh+amber_data->
Nbona; ++i)
10718 { bonds[numBonds].atom1 = amber_data->
BondAt1[i-amber_data->
Nbonh] / 3;
10719 bonds[numBonds].atom2 = amber_data->
BondAt2[i-amber_data->
Nbonh] / 3;
10720 bonds[numBonds].bond_type = amber_data->
BondNum[i-amber_data->
Nbonh] - 1;
10721 if (bonds[i].atom1>=numAtoms || bonds[i].atom2>=numAtoms ||
10722 bonds[i].bond_type>=amber_data->
Numbnd)
10723 {
char err_msg[128];
10724 sprintf(err_msg,
"BOND (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-amber_data->
Nbonh);
10733 if ( numBonds != amber_data->
Nbonh + amber_data->
Nbona) {
10735 " bonds with zero force constants.\n" <<
endi;
10737 "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
endi;
10743 { ntheth = amber_data->
Ntheth;
10744 angles =
new Angle[numAngles];
10745 if (angles == NULL || numAngles < ntheth)
10746 NAMD_die(
"memory allocation failed when reading angle information");
10748 for (i=0; i<ntheth; ++i)
10749 { angles[i].atom1 = amber_data->
AngleHAt1[i] / 3;
10750 angles[i].atom2 = amber_data->
AngleHAt2[i] / 3;
10751 angles[i].atom3 = amber_data->
AngleHAt3[i] / 3;
10752 angles[i].angle_type = amber_data->
AngleHNum[i] - 1;
10753 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
10754 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
10755 {
char err_msg[128];
10756 sprintf(err_msg,
"ANGLE (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
10761 for (i=ntheth; i<numAngles; ++i)
10762 { angles[i].atom1 = amber_data->
AngleAt1[i-ntheth] / 3;
10763 angles[i].atom2 = amber_data->
AngleAt2[i-ntheth] / 3;
10764 angles[i].atom3 = amber_data->
AngleAt3[i-ntheth] / 3;
10765 angles[i].angle_type = amber_data->
AngleNum[i-ntheth] - 1;
10766 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
10767 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
10768 {
char err_msg[128];
10769 sprintf(err_msg,
"ANGLE (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-ntheth);
10784 if (exclusions == NULL && amber_data->
Nnb > 0)
10785 NAMD_die(
"memory allocation failed when reading exclusion information");
10787 for (i=0; i<numAtoms; ++i)
10788 for (j=0; j<amber_data->
Iblo[i]; ++j)
10789 {
if (current_index >= amber_data->
Nnb)
10790 {
char err_msg[128];
10791 sprintf(err_msg,
"EXCLUSION INDEX EXCEEDS NUMBER OF EXLCUSIONS %d IN AMBER FILE, AT ATOM #%d\n",
10792 amber_data->
Nnb, i+1);
10797 if (amber_data->
ExclAt[current_index] != 0)
10800 a2 = amber_data->
ExclAt[current_index] - 1;
10806 sprintf(err_msg,
"Atom #%d has exclusion with atom #%d, in reverse order.", i+1, a2+1);
10810 {
char err_msg[128];
10811 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN AMBER FILE\n", i+1);
10814 else if (a2 >= numAtoms)
10815 {
char err_msg[128];
10816 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN AMBER FILE",
10817 a2+1, numAtoms, current_index+1);
10820 exclusions[numExclusions].atom1 = i;
10821 exclusions[numExclusions].atom2 = a2;
10826 if (current_index < amber_data->Nnb)
10827 {
char err_msg[128];
10828 sprintf(err_msg,
"Num of exclusions recorded (%d) is smaller than what it's supposed to be (%d)",
10829 current_index,amber_data->
Nnb);
10835 numDihedrals = amber_data->
Nphih + amber_data->
Nphia;
10836 if (numDihedrals > 0)
10837 { nphih = amber_data->
Nphih;
10838 dihedrals =
new Dihedral[numDihedrals];
10839 if (dihedrals == NULL || numDihedrals < nphih)
10840 NAMD_die(
"memory allocation failed when reading dihedral information");
10842 for (i=0; i<nphih; ++i)
10843 { dihedrals[i].atom1 = amber_data->
DihHAt1[i] / 3;
10844 dihedrals[i].atom2 = amber_data->
DihHAt2[i] / 3;
10845 dihedrals[i].atom3 = amber_data->
DihHAt3[i] / 3;
10846 dihedrals[i].atom4 = amber_data->
DihHAt4[i] / 3;
10847 dihedrals[i].dihedral_type = amber_data->
DihHNum[i] - 1;
10850 for (i=nphih; i<numDihedrals; ++i)
10851 { dihedrals[i].atom1 = amber_data->
DihAt1[i-nphih] / 3;
10852 dihedrals[i].atom2 = amber_data->
DihAt2[i-nphih] / 3;
10853 dihedrals[i].atom3 = amber_data->
DihAt3[i-nphih] / 3;
10854 dihedrals[i].atom4 = amber_data->
DihAt4[i-nphih] / 3;
10855 dihedrals[i].dihedral_type = amber_data->
DihNum[i-nphih] - 1;
10870 for (i=0; i<numDihedrals; ++i)
10871 {
if (dihedrals[i].atom3 < 0 || dihedrals[i].atom4 < 0)
10872 { dihedrals[i].atom3 = abs(dihedrals[i].atom3);
10873 dihedrals[i].atom4 = abs(dihedrals[i].atom4);
10876 {
if (dihedrals[i].atom1 < dihedrals[i].atom4)
10877 a1=dihedrals[i].atom1, a2=dihedrals[i].atom4;
10879 a1=dihedrals[i].atom4, a2=dihedrals[i].atom1;
10883 min=0, max=numExclusions-1;
10884 while (!found && min<=max)
10885 { index = (min+max)/2;
10886 if (exclusions[index].atom1 == a1)
10888 else if (exclusions[index].atom1 < a1)
10894 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
10897 for (j=index-1; j>=0 && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; --j);
10898 if (j<0 || exclusions[j].atom1!=a1)
10899 for (j=index; j<numExclusions && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; ++j);
10900 if (j<numExclusions && exclusions[j].atom1==a1)
10901 exclusions[j].modified = 1;
10903 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
10905 if (dihedrals[i].atom1>=numAtoms || dihedrals[i].atom2>=numAtoms ||
10906 dihedrals[i].atom3>=numAtoms || dihedrals[i].atom4>=numAtoms ||
10907 dihedrals[i].dihedral_type>=amber_data->
Nptra)
10908 {
char err_msg[128];
10909 sprintf(err_msg,
"DIHEDRAL # %d OVERFLOW IN PARM FILE", i+1);
10915 if (amber_data -> HasCMAP) {
10917 if (numCrossterms > 0) {
10918 crossterms =
new Crossterm[numCrossterms];
10920 for (i = 0; i < numCrossterms; ++i) {
10923 crossterms[i].atom2 = amber_data->
CMAPIndex[6*i+1]-1;
10924 crossterms[i].atom3 = amber_data->
CMAPIndex[6*i+2]-1;
10925 crossterms[i].atom4 = amber_data->
CMAPIndex[6*i+3]-1;
10927 crossterms[i].atom5 = amber_data->
CMAPIndex[6*i+1]-1;
10928 crossterms[i].atom6 = amber_data->
CMAPIndex[6*i+2]-1;
10929 crossterms[i].atom7 = amber_data->
CMAPIndex[6*i+3]-1;
10930 crossterms[i].atom8 = amber_data->
CMAPIndex[6*i+4]-1;
10932 crossterms[i].crossterm_type = amber_data->
CMAPIndex[6*i+5]-1;
10937 numRealBonds = numBonds;
10938 build_atom_status();
10954 void Molecule::read_parm(
Ambertoppar *amber_data)
10956 #ifdef MEM_OPT_VERSION
10957 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
10959 int i,j,ntheth,nphih,current_index,a1,a2,
10960 max,min,index,found;
10963 NAMD_die(
"No data read from parm file yet!");
10966 numAtoms = amber_data->
Natom;
10967 atoms =
new Atom[numAtoms];
10974 if (atoms == NULL || atomNames == NULL )
10975 NAMD_die(
"memory allocation failed when reading atom information");
10977 for (i=0; i<numAtoms; ++i)
10978 { atomNames[i].resname = nameArena->getNewArray(5);
10979 atomNames[i].atomname = nameArena->getNewArray(5);
10980 atomNames[i].atomtype = nameArena->getNewArray(5);
10981 if (atomNames[i].resname == NULL || atomNames[i].atomname == NULL || atomNames[i].atomtype == NULL)
10982 NAMD_die(
"memory allocation failed when reading atom information");
10983 for (j=0; j<4; ++j)
10984 { atomNames[i].resname[j] = amber_data->
ResNames[amber_data->
AtomRes[i]*4+j];
10985 atomNames[i].atomname[j] = amber_data->
AtomNames[i*4+j];
10986 atomNames[i].atomtype[j] = amber_data->
AtomSym[i*4+j];
10988 atomNames[i].resname[4] = atomNames[i].atomname[4] = atomNames[i].atomtype[4] =
'\0';
10989 strtok(atomNames[i].resname,
" ");
10990 strtok(atomNames[i].atomname,
" ");
10991 strtok(atomNames[i].atomtype,
" ");
10992 atoms[i].mass = amber_data->
Masses[i];
10994 atoms[i].charge = amber_data->
Charges[i] / 18.2223;
10995 atoms[i].vdw_type = amber_data->
Iac[i] - 1;
10998 if ( tmpResLookup ) tmpResLookup =
11001 if(atomSegResids) {
11003 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
11011 }
else if (atoms[i].mass <= 0.05) {
11012 ++numZeroMassAtoms;
11014 }
else if (atoms[i].mass < 1.0) {
11016 }
else if (atoms[i].mass <=3.5) {
11018 }
else if ((atomNames[i].atomname[0] ==
'O') &&
11019 (atoms[i].mass >= 14.0) &&
11020 (atoms[i].mass <= 18.0)) {
11033 if (amber_data->
Nbonh + amber_data->
Nbona > 0)
11035 if (bonds == NULL || amber_data->
Nbona < 0)
11036 NAMD_die(
"memory allocation failed when reading bond information");
11038 for (i=0; i<amber_data->
Nbonh; ++i)
11039 { bonds[numBonds].atom1 = amber_data->
BondHAt1[i] / 3;
11040 bonds[numBonds].atom2 = amber_data->
BondHAt2[i] / 3;
11041 bonds[numBonds].bond_type = amber_data->
BondHNum[i] - 1;
11042 if (bonds[numBonds].atom1>=numAtoms || bonds[numBonds].atom2>=numAtoms ||
11043 bonds[numBonds].bond_type>=amber_data->
Numbnd)
11044 {
char err_msg[128];
11045 sprintf(err_msg,
"BOND (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11053 for (i=amber_data->
Nbonh; i<amber_data->Nbonh+amber_data->
Nbona; ++i)
11054 { bonds[numBonds].atom1 = amber_data->
BondAt1[i-amber_data->
Nbonh] / 3;
11055 bonds[numBonds].atom2 = amber_data->
BondAt2[i-amber_data->
Nbonh] / 3;
11056 bonds[numBonds].bond_type = amber_data->
BondNum[i-amber_data->
Nbonh] - 1;
11057 if (bonds[i].atom1>=numAtoms || bonds[i].atom2>=numAtoms ||
11058 bonds[i].bond_type>=amber_data->
Numbnd)
11059 {
char err_msg[128];
11060 sprintf(err_msg,
"BOND (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-amber_data->
Nbonh);
11069 if ( numBonds != amber_data->
Nbonh + amber_data->
Nbona) {
11071 " bonds with zero force constants.\n" <<
endi;
11073 "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
endi;
11079 { ntheth = amber_data->
Ntheth;
11080 angles =
new Angle[numAngles];
11081 if (angles == NULL || numAngles < ntheth)
11082 NAMD_die(
"memory allocation failed when reading angle information");
11084 for (i=0; i<ntheth; ++i)
11085 { angles[i].atom1 = amber_data->
AngleHAt1[i] / 3;
11086 angles[i].atom2 = amber_data->
AngleHAt2[i] / 3;
11087 angles[i].atom3 = amber_data->
AngleHAt3[i] / 3;
11088 angles[i].angle_type = amber_data->
AngleHNum[i] - 1;
11089 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11090 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11091 {
char err_msg[128];
11092 sprintf(err_msg,
"ANGLE (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11097 for (i=ntheth; i<numAngles; ++i)
11098 { angles[i].atom1 = amber_data->
AngleAt1[i-ntheth] / 3;
11099 angles[i].atom2 = amber_data->
AngleAt2[i-ntheth] / 3;
11100 angles[i].atom3 = amber_data->
AngleAt3[i-ntheth] / 3;
11101 angles[i].angle_type = amber_data->
AngleNum[i-ntheth] - 1;
11102 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11103 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11104 {
char err_msg[128];
11105 sprintf(err_msg,
"ANGLE (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-ntheth);
11120 if (exclusions == NULL && amber_data->
Nnb > 0)
11121 NAMD_die(
"memory allocation failed when reading exclusion information");
11123 for (i=0; i<numAtoms; ++i)
11124 for (j=0; j<amber_data->
Iblo[i]; ++j)
11125 {
if (current_index >= amber_data->
Nnb)
11126 {
char err_msg[128];
11127 sprintf(err_msg,
"EXCLUSION INDEX EXCEEDS NUMBER OF EXLCUSIONS %d IN AMBER FILE, AT ATOM #%d\n",
11128 amber_data->
Nnb, i+1);
11133 if (amber_data->
ExclAt[current_index] != 0)
11136 a2 = amber_data->
ExclAt[current_index] - 1;
11142 sprintf(err_msg,
"Atom #%d has exclusion with atom #%d, in reverse order.", i+1, a2+1);
11146 {
char err_msg[128];
11147 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN AMBER FILE\n", i+1);
11150 else if (a2 >= numAtoms)
11151 {
char err_msg[128];
11152 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN AMBER FILE",
11153 a2+1, numAtoms, current_index+1);
11156 exclusions[numExclusions].atom1 = i;
11157 exclusions[numExclusions].atom2 = a2;
11162 if (current_index < amber_data->Nnb)
11163 {
char err_msg[128];
11164 sprintf(err_msg,
"Num of exclusions recorded (%d) is smaller than what it's supposed to be (%d)",
11165 current_index,amber_data->
Nnb);
11171 numDihedrals = amber_data->
Nphih + amber_data->
Nphia;
11172 if (numDihedrals > 0)
11173 { nphih = amber_data->
Nphih;
11174 dihedrals =
new Dihedral[numDihedrals];
11175 if (dihedrals == NULL || numDihedrals < nphih)
11176 NAMD_die(
"memory allocation failed when reading dihedral information");
11178 for (i=0; i<nphih; ++i)
11179 { dihedrals[i].atom1 = amber_data->
DihHAt1[i] / 3;
11180 dihedrals[i].atom2 = amber_data->
DihHAt2[i] / 3;
11181 dihedrals[i].atom3 = amber_data->
DihHAt3[i] / 3;
11182 dihedrals[i].atom4 = amber_data->
DihHAt4[i] / 3;
11183 dihedrals[i].dihedral_type = amber_data->
DihHNum[i] - 1;
11186 for (i=nphih; i<numDihedrals; ++i)
11187 { dihedrals[i].atom1 = amber_data->
DihAt1[i-nphih] / 3;
11188 dihedrals[i].atom2 = amber_data->
DihAt2[i-nphih] / 3;
11189 dihedrals[i].atom3 = amber_data->
DihAt3[i-nphih] / 3;
11190 dihedrals[i].atom4 = amber_data->
DihAt4[i-nphih] / 3;
11191 dihedrals[i].dihedral_type = amber_data->
DihNum[i-nphih] - 1;
11206 for (i=0; i<numDihedrals; ++i)
11207 {
if (dihedrals[i].atom3 < 0 || dihedrals[i].atom4 < 0)
11208 { dihedrals[i].atom3 = abs(dihedrals[i].atom3);
11209 dihedrals[i].atom4 = abs(dihedrals[i].atom4);
11212 {
if (dihedrals[i].atom1 < dihedrals[i].atom4)
11213 a1=dihedrals[i].atom1, a2=dihedrals[i].atom4;
11215 a1=dihedrals[i].atom4, a2=dihedrals[i].atom1;
11219 min=0, max=numExclusions-1;
11220 while (!found && min<=max)
11221 { index = (min+max)/2;
11222 if (exclusions[index].atom1 == a1)
11224 else if (exclusions[index].atom1 < a1)
11230 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11233 for (j=index-1; j>=0 && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; --j);
11234 if (j<0 || exclusions[j].atom1!=a1)
11235 for (j=index; j<numExclusions && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; ++j);
11236 if (j<numExclusions && exclusions[j].atom1==a1)
11237 exclusions[j].modified = 1;
11239 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11241 if (dihedrals[i].atom1>=numAtoms || dihedrals[i].atom2>=numAtoms ||
11242 dihedrals[i].atom3>=numAtoms || dihedrals[i].atom4>=numAtoms ||
11243 dihedrals[i].dihedral_type>=amber_data->
Nptra)
11244 {
char err_msg[128];
11245 sprintf(err_msg,
"DIHEDRAL # %d OVERFLOW IN PARM FILE", i+1);
11251 numRealBonds = numBonds;
11252 build_atom_status();
11269 initialize(simParams,param);
11271 read_parm(gromacsTopFile);
11273 #ifndef MEM_OPT_VERSION
11276 assignLCPOTypes( 3 );
11294 #ifdef MEM_OPT_VERSION
11295 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
11303 atoms =
new Atom[numAtoms];
11310 if (atoms == NULL || atomNames == NULL )
11311 NAMD_die(
"memory allocation failed when reading atom information");
11315 for (i=0; i<numAtoms; ++i) {
11316 char *resname = nameArena->getNewArray(11);
11317 char *atomname = nameArena->getNewArray(11);
11318 char *atomtype = nameArena->getNewArray(11);
11319 int resnum,typenum;
11322 if (resname == NULL || atomname == NULL || atomtype == NULL)
11323 NAMD_die(
"memory allocation failed when reading atom information");
11326 gf->
getAtom(i,&resnum,resname,atomname,atomtype,&typenum,
11329 atomNames[i].resname = resname;
11330 atomNames[i].atomname = atomname;
11331 atomNames[i].atomtype = atomtype;
11332 atoms[i].mass = mass;
11333 atoms[i].charge =
charge;
11334 atoms[i].vdw_type = typenum;
11337 if ( tmpResLookup ) tmpResLookup =
11338 tmpResLookup->
append(
"MAIN", resnum+1, i);
11340 if(atomSegResids) {
11342 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
11343 one->
resid = resnum+1;
11352 }
else if (atoms[i].mass <= 0.05) {
11354 }
else if (atoms[i].mass < 1.0) {
11356 }
else if (atoms[i].mass <=3.5) {
11358 }
else if ((atomNames[i].atomname[0] ==
'O') &&
11359 (atoms[i].mass >= 14.0) &&
11360 (atoms[i].mass <= 18.0)) {
11367 bonds =
new Bond[numBonds];
11369 NAMD_die(
"memory allocation failed when reading bond information");
11370 for(i=0;i<numBonds;i++) {
11373 gf->
getBond(i,&atom1,&atom2,&type);
11374 bonds[i].atom1 = atom1;
11375 bonds[i].atom2 = atom2;
11376 bonds[i].bond_type = (
Index)type;
11381 angles =
new Angle[numAngles];
11382 if (angles == NULL)
11383 NAMD_die(
"memory allocation failed when reading angle information");
11384 for(i=0;i<numAngles;i++) {
11386 int atom1,atom2,atom3;
11387 gf->
getAngle(i,&atom1,&atom2,&atom3,&type);
11389 angles[i].atom1 = atom1;
11390 angles[i].atom2 = atom2;
11391 angles[i].atom3 = atom3;
11393 angles[i].angle_type=type;
11397 exclusions =
new Exclusion[numExclusions];
11461 dihedrals =
new Dihedral[numDihedrals];
11462 if (dihedrals == NULL)
11463 NAMD_die(
"memory allocation failed when reading dihedral information");
11464 for(i=0;i<numDihedrals;i++) {
11466 int atom1,atom2,atom3,atom4;
11467 gf->
getDihedral(i,&atom1,&atom2,&atom3,&atom4,&type);
11468 dihedrals[i].atom1 = atom1;
11469 dihedrals[i].atom2 = atom2;
11470 dihedrals[i].atom3 = atom3;
11471 dihedrals[i].atom4 = atom4;
11472 dihedrals[i].dihedral_type = type;
11485 const_cast<GromacsTopFile*
>(gf)->getPairLJArrays2(indxLJA, indxLJB, pairC6, pairC12);
11488 gromacsPair_type[i] = i;
11497 pointerToLJBeg =
new int[numAtoms];
11498 pointerToLJEnd =
new int[numAtoms];
11500 for(
int i=0; i < numAtoms; i++) {
11501 pointerToLJBeg[i] = -1;
11502 pointerToLJEnd[i] = -2;
11505 if(pointerToLJBeg[indxLJA[i]] == -1) {
11506 pointerToLJBeg[indxLJA[i]] = i;
11507 oldIndex = indxLJA[i];
11509 pointerToLJEnd[oldIndex] = i;
11522 const_cast<GromacsTopFile*
>(gf)->getPairGaussArrays2(indxGaussA, indxGaussB, gA, gMu1, giSigma1, gMu2, giSigma2, gRepulsive);
11525 pointerToGaussBeg =
new int[numAtoms];
11526 pointerToGaussEnd =
new int[numAtoms];
11527 for(
int i=0; i < numAtoms; i++) {
11528 pointerToGaussBeg[i] = -1;
11529 pointerToGaussEnd[i] = -2;
11533 if(pointerToGaussBeg[indxGaussA[i]] == -1) {
11534 pointerToGaussBeg[indxGaussA[i]] = i;
11535 oldIndex = indxGaussA[i];
11537 pointerToGaussEnd[oldIndex] = i;
11540 iout <<
iINFO <<
"Finished reading explicit pair from Gromacs file:\n" <<
11541 iINFO <<
"Found a total of: " << numPair <<
" explicit pairs--of which: " <<
11542 numLJPair <<
" are LJ style pairs and " << numGaussPair <<
11543 " are Gaussian style pairs.\n" <<
endi;
11547 #if GROMACS_EXCLUSIONS
11550 int* atom1 =
new int[numExclusions];
11551 int* atom2 =
new int[numExclusions];
11552 for(
int j=0; j<numExclusions;j++) {
11558 read_exclusions(atom1,atom2,numExclusions);
11622 numRealBonds = numBonds;
11623 build_atom_status();
11628 #ifndef MEM_OPT_VERSION
11642 #ifdef MEM_OPT_VERSION
11644 Index Molecule::insert_new_mass(
Real newMass){
11646 for(
int i=massPoolSize-1; i>=0; i--){
11647 if(fabs(atomMassPool[i]-newMass)<=1e-6)
11651 Real *tmp =
new Real[massPoolSize+1];
11652 tmp[massPoolSize] = newMass;
11653 memcpy((
void *)tmp, (
const void *)atomMassPool,
sizeof(
Real)*massPoolSize);
11654 delete [] atomMassPool;
11655 atomMassPool = tmp;
11657 return (
Index)(massPoolSize-1);
11660 void Molecule::addNewExclSigPool(
const vector<ExclusionSignature>& newExclSigPool){
11662 for(
int i=0; i<exclSigPoolSize; i++)
11663 tmpExclSigPool[i] = exclSigPool[i];
11664 for(
int i=0; i<newExclSigPool.size(); i++)
11665 tmpExclSigPool[i+exclSigPoolSize] = newExclSigPool[i];
11667 exclSigPoolSize += newExclSigPool.size();
11668 exclSigPool = tmpExclSigPool;
11672 msg->
put((
short)tupleType);
11673 msg->
put(numOffset);
11674 msg->
put(numOffset, offset);
11675 msg->
put(tupleParamType);
11684 msg->
get(numOffset);
11686 offset =
new int[numOffset];
11687 msg->
get(numOffset*
sizeof(
int), (
char *)offset);
11689 msg->
get(tupleParamType);
11695 for(
int i=0; i<bondCnt; i++)
11696 bondSigs[i].pack(msg);
11698 msg->
put(angleCnt);
11699 for(
int i=0; i<angleCnt; i++)
11700 angleSigs[i].pack(msg);
11702 msg->
put(dihedralCnt);
11703 for(
int i=0; i<dihedralCnt; i++)
11704 dihedralSigs[i].pack(msg);
11706 msg->
put(improperCnt);
11707 for(
int i=0; i<improperCnt; i++)
11708 improperSigs[i].pack(msg);
11710 msg->
put(crosstermCnt);
11711 for(
int i=0; i<crosstermCnt; i++)
11712 crosstermSigs[i].pack(msg);
11715 msg->
put(gromacsPairCnt);
11716 for(
int i=0; i<gromacsPairCnt; i++)
11717 gromacsPairSigs[i].pack(msg);
11722 delete [] bondSigs;
11725 for(
int i=0; i<bondCnt; i++)
11726 bondSigs[i].unpack(msg);
11727 }
else bondSigs = NULL;
11729 msg->
get(angleCnt);
11730 delete [] angleSigs;
11733 for(
int i=0; i<angleCnt; i++)
11734 angleSigs[i].unpack(msg);
11735 }
else angleSigs = NULL;
11737 msg->
get(dihedralCnt);
11738 delete [] dihedralSigs;
11741 for(
int i=0; i<dihedralCnt; i++)
11742 dihedralSigs[i].unpack(msg);
11743 }
else dihedralSigs = NULL;
11745 msg->
get(improperCnt);
11746 delete [] improperSigs;
11749 for(
int i=0; i<improperCnt; i++)
11750 improperSigs[i].unpack(msg);
11751 }
else improperSigs = NULL;
11753 msg->
get(crosstermCnt);
11754 delete [] crosstermSigs;
11755 if(crosstermCnt>0){
11757 for(
int i=0; i<crosstermCnt; i++)
11758 crosstermSigs[i].unpack(msg);
11759 }
else crosstermSigs = NULL;
11763 msg->
get(gromacsPairCnt);
11764 delete [] gromacsPairSigs;
11765 if(gromacsPairCnt>0){
11767 for(
int i=0; i<gromacsPairCnt; i++)
11768 gromacsPairSigs[i].unpack(msg);
11769 }
else gromacsPairSigs = NULL;
11783 origTupleCnt = bondCnt;
11784 tupleSigs= bondSigs;
11785 for(
int i=0; i<origTupleCnt; i++){
11786 if(tupleSigs[i].isEmpty())
11790 delete [] tupleSigs;
11792 }
else if(bondCnt!=origTupleCnt){
11795 for(
int i=0; i<origTupleCnt; i++){
11796 if(!tupleSigs[i].isEmpty()){
11797 newTupleSigs[idx] = tupleSigs[i];
11801 delete [] tupleSigs;
11802 bondSigs = newTupleSigs;
11808 origTupleCnt = angleCnt;
11809 tupleSigs = angleSigs;
11810 for(
int i=0; i<origTupleCnt; i++){
11811 if(tupleSigs[i].isEmpty())
11815 delete [] tupleSigs;
11817 }
else if(angleCnt!=origTupleCnt){
11820 for(
int i=0; i<origTupleCnt; i++){
11821 if(!tupleSigs[i].isEmpty()){
11822 newTupleSigs[idx] = tupleSigs[i];
11826 delete [] tupleSigs;
11827 angleSigs = newTupleSigs;
11833 origTupleCnt = dihedralCnt;
11834 tupleSigs = dihedralSigs;
11835 for(
int i=0; i<origTupleCnt; i++){
11836 if(tupleSigs[i].isEmpty())
11839 if(dihedralCnt==0){
11840 delete [] tupleSigs;
11841 dihedralSigs = NULL;
11842 }
else if(dihedralCnt!=origTupleCnt){
11845 for(
int i=0; i<origTupleCnt; i++){
11846 if(!tupleSigs[i].isEmpty()){
11847 newTupleSigs[idx] = tupleSigs[i];
11851 delete [] tupleSigs;
11852 dihedralSigs = newTupleSigs;
11859 origTupleCnt = improperCnt;
11860 tupleSigs = improperSigs;
11861 for(
int i=0; i<origTupleCnt; i++){
11862 if(tupleSigs[i].isEmpty())
11865 if(improperCnt==0){
11866 delete [] tupleSigs;
11867 improperSigs = NULL;
11868 }
else if(improperCnt!=origTupleCnt){
11871 for(
int i=0; i<origTupleCnt; i++){
11872 if(!tupleSigs[i].isEmpty()){
11873 newTupleSigs[idx] = tupleSigs[i];
11877 delete [] tupleSigs;
11878 improperSigs = newTupleSigs;
11884 origTupleCnt = crosstermCnt;
11885 tupleSigs = crosstermSigs;
11886 for(
int i=0; i<origTupleCnt; i++){
11887 if(tupleSigs[i].isEmpty())
11890 if(crosstermCnt==0){
11891 delete [] tupleSigs;
11892 crosstermSigs = NULL;
11893 }
else if(crosstermCnt!=origTupleCnt){
11896 for(
int i=0; i<origTupleCnt; i++){
11897 if(!tupleSigs[i].isEmpty()){
11898 newTupleSigs[idx] = tupleSigs[i];
11902 delete [] tupleSigs;
11903 crosstermSigs = newTupleSigs;
11910 origTupleCnt = gromacsPairCnt;
11911 tupleSigs = gromacsPairSigs;
11912 for(
int i=0; i<origTupleCnt; i++){
11913 if(tupleSigs[i].isEmpty())
11916 if(gromacsPairCnt==0){
11917 delete [] tupleSigs;
11918 gromacsPairSigs = NULL;
11919 }
else if(gromacsPairCnt!=origTupleCnt){
11922 for(
int i=0; i<origTupleCnt; i++){
11923 if(!tupleSigs[i].isEmpty()){
11924 newTupleSigs[idx] = tupleSigs[i];
11928 delete [] tupleSigs;
11929 gromacsPairSigs = newTupleSigs;
11939 for(
int i=0; i<fullExclCnt; i++){
11940 if(fullOffset[i]==0)
continue;
11945 delete [] fullOffset;
11947 }
else if(newCnt!=fullExclCnt){
11948 int *tmpOffset =
new int[newCnt];
11950 for(
int i=0; i<fullExclCnt; i++){
11951 if(fullOffset[i]==0)
continue;
11952 tmpOffset[newCnt] = fullOffset[i];
11955 delete [] fullOffset;
11956 fullOffset = tmpOffset;
11957 fullExclCnt = newCnt;
11962 for(
int i=0; i<modExclCnt; i++){
11963 if(modOffset[i]==0)
continue;
11968 delete [] modOffset;
11970 }
else if(newCnt!=modExclCnt){
11971 int *tmpOffset =
new int[newCnt];
11973 for(
int i=0; i<modExclCnt; i++){
11974 if(modOffset[i]==0)
continue;
11975 tmpOffset[newCnt] = modOffset[i];
11978 delete [] modOffset;
11979 modOffset = tmpOffset;
11980 modExclCnt = newCnt;
11995 int high = fullExclCnt-1;
11996 int mid = (low+high)/2;
11998 if(offset<fullOffset[mid]){
12000 mid = (high+low)/2;
12001 }
else if(offset>fullOffset[mid]){
12003 mid = (high+low)/2;
12009 if(retidx!=-1)
return retidx;
12013 high = modExclCnt-1;
12014 mid = (low+high)/2;
12016 if(offset<modOffset[mid]){
12018 mid = (high+low)/2;
12019 }
else if(offset>modOffset[mid]){
12021 mid = (high+low)/2;
12031 msg->
put(fullExclCnt);
12032 msg->
put(fullExclCnt, fullOffset);
12033 msg->
put(modExclCnt);
12034 msg->
put(modExclCnt, modOffset);
12038 msg->
get(fullExclCnt);
12039 delete [] fullOffset;
12040 fullOffset =
new int[fullExclCnt];
12041 msg->
get(fullExclCnt*
sizeof(
int), (
char *)fullOffset);
12042 msg->
get(modExclCnt);
12043 delete [] modOffset;
12044 modOffset =
new int[modExclCnt];
12045 msg->
get(modExclCnt*
sizeof(
int), (
char *)modOffset);
12046 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
12052 #endif // MOLECULE2_C defined = second object file
int get_atom_from_name(const char *segid, int resid, const char *aname) const
int get_residue_size(const char *segid, int resid) const
void build_gridforce_params(StringList *, StringList *, StringList *, StringList *, PDB *, char *)
void getBond(int num, int *atomi, int *atomj, int *bondtype) const
void assign_improper_index(const char *, const char *, const char *, const char *, Improper *, int)
std::ostream & iINFO(std::ostream &s)
void flipNum(char *elem, int elemSize, int numElems)
#define COMPRESSED_PSF_MAGICNUM
void NAMD_err(const char *err_msg)
void print_bonds(Parameters *)
int NAMD_read_line(FILE *fd, char *buf, int bufsize)
void assign_vdw_index(const char *, Atom *)
unsigned int hydrogenGroupSize
void build_extra_bonds(Parameters *parameters, StringList *file)
void setOccupancyData(molfile_atom_t *atomarray)
void read_alch_unpert_dihedrals(FILE *)
void assign_bond_index(const char *, const char *, Bond *)
vector< string > ResNames
Bool is_hydrogenGroupParent(int)
static GridforceGrid * new_grid(int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
void send_Molecule(MOStream *)
static __thread unsigned int * exclusions
static __thread atom * atoms
void unpack(MIStream *msg)
static void pack_grid(GridforceGrid *grid, MOStream *msg)
std::ostream & endi(std::ostream &s)
static __thread int2 * exclusionsByAtom
int get_atom_from_index_in_residue(const char *segid, int resid, int index) const
std::ostream & iWARN(std::ostream &s)
void setBFactorData(molfile_atom_t *atomarray)
MIStream * get(char &data)
DihedralValue * dihedral_array
int add(const Elem &elem)
void reloadCharges(float charge[], int n)
#define NAMD_FILENAME_BUFFER_SIZE
MGridforceParams * get_first()
int getLCPOTypeAmber(char atomType[11], int numBonds)
void unpack(MIStream *msg)
int getNumExclusions() const
#define COMPRESSED_PSF_VER
FourBodyConsts values[MAX_MULTIPLICITY]
void build_langevin_params(BigReal coupling, BigReal drudeCoupling, Bool doHydrogen)
void getDihedral(int num, int *atomi, int *atomj, int *atomk, int *atoml, int *type) const
static Units next(Units u)
void get_angle_params(Real *k, Real *theta0, Real *k_ub, Real *r_ub, Index index)
BigReal getEnergyTailCorr(const BigReal, const int)
void assign_crossterm_index(const char *, const char *, const char *, const char *, const char *, const char *, const char *, const char *, Crossterm *)
PDBAtom * atom(int place)
void NAMD_bug(const char *err_msg)
int NAMD_blank_string(char *str)
#define SPLIT_PATCH_HYDROGEN
int getNumGaussPair() const
int checkexcl(int atom1, int atom2) const
void build_constraint_params(StringList *, StringList *, StringList *, PDB *, char *)
int findOffset(int offset, int *fullOrMod)
void build_constorque_params(StringList *, StringList *, StringList *, StringList *, StringList *, StringList *, PDB *, char *)
void read_alch_unpert_bonds(FILE *)
FILE * Fopen(const char *filename, const char *mode)
int lookup(const char *segid, int resid, int *begin, int *end) const
int NAMD_find_word(const char *source, const char *search)
void build_stirred_atoms(StringList *, StringList *, PDB *, char *)
ImproperValue * improper_array
int get_vdw_pair_params(Index ind1, Index ind2, Real *, Real *, Real *, Real *)
void removeEmptyTupleSigs()
void NAMD_die(const char *err_msg)
HashPool< HashString > atomNamePool
BigReal getVirialTailCorr(const BigReal)
void build_ss_flags(const StringList *ssfile, const StringList *sscol, PDB *initial_pdb, const char *cwd)
void build_rotdrag_params(StringList *, StringList *, StringList *, StringList *, StringList *, StringList *, PDB *, char *)
HashPool< HashString > resNamePool
int get_num_vdw_params(void)
FourBodyConsts values[MAX_MULTIPLICITY]
MGridforceParamsList mgridforcelist
void build_fep_flags(StringList *, StringList *, PDB *, char *, const char *)
vector< string > AtomNames
ResidueLookupElem * append(const char *segid, int resid, int aid)
int add(const Elem &elem)
void build_alch_unpert_bond_lists(char *)
struct OutputAtomRecord::floatVals fSet
void build_fixed_atoms(StringList *, StringList *, PDB *, char *)
BlockRadixSort::TempStorage sort
void compute_LJcorrection()
int pressureProfileAtomTypes
int atomsInMigrationGroup
void print_atoms(Parameters *)
static GridforceGrid * unpack_grid(int gridnum, MIStream *msg)
struct OutputAtomRecord::integerVals iSet
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ cudaTextureObject_t const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t float const PatchPairRecord *__restrict__ const int *__restrict__ atomIndex
void getAngle(int num, int *atomi, int *atomj, int *atomk, int *angletype) const
Bool extraBondsCosAnglesSetByUser
int atomsInMigrationGroup
void assign_dihedral_index(const char *, const char *, const char *, const char *, Dihedral *, int, int)
int NAMD_read_int(FILE *fd, const char *msg)
HashPool< HashString > segNamePool
void delete_alch_bonded(void)
StringList * find(const char *name) const
HashPool< HashString > atomTypePool
void build_exPressure_atoms(StringList *, StringList *, PDB *, char *)
void build_constant_forces(char *)
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
void unpack(MIStream *msg)
void assign_angle_index(const char *, const char *, const char *, Angle *, int)
BigReal getVdwLambda(const BigReal)
MOStream * put(char data)
int pairInteractionGroup2
int getNumDihedrals() const
std::ostream & iERROR(std::ostream &s)
void get_bond_params(Real *k, Real *x0, Index index)
int pairInteractionGroup1
void receive_Molecule(MIStream *)
void build_movdrag_params(StringList *, StringList *, StringList *, PDB *, char *)
ExclusionSettings exclude
struct OutputAtomRecord::shortVals sSet
void get_vdw_params(Real *sigma, Real *epsilon, Real *sigma14, Real *epsilon14, Index index)
int getLCPOTypeCharmm(char atomType[11], int numBonds)
ResizeArray< int > atomIndex
void getAtom(int num, int *residue_number, char *residue_name, char *atom_name, char *atom_type, int *atom_typenum, Real *charge, Real *mass) const
Molecule(SimParameters *, Parameters *param)
HashPool< AtomSigInfo > atomSigPool
void read_alch_unpert_angles(FILE *)