49 #define MIN_DEBUG_LEVEL 3 61 #define M_PI 3.14159265358979323846 65 #define GROMACS_PAIR 1 68 #ifndef GROMACS_EXCLUSIONS 69 #define GROMACS_EXCLUSIONS 1 74 #ifndef MOLECULE2_C // first object file 76 #ifdef MEM_OPT_VERSION 77 template int lookupCstPool<AtomSignature>(
const vector<AtomSignature>&,
const AtomSignature&);
78 template int lookupCstPool<ExclusionSignature>(
const vector<ExclusionSignature>&,
const ExclusionSignature&);
82 const char *segid,
int resid,
int *begin,
int *end)
const {
85 while ( elem && strcasecmp(elem->
mySegid,segid) ) elem = elem->
next;
86 if ( elem && (resid >= elem->
firstResid) && (resid <= elem->lastResid) ) {
95 const char *segid,
int resid,
int aid) {
97 if ( firstResid == -1 ) {
98 strcpy(mySegid,segid);
102 atomIndex.add(aid+1);
103 }
else if ( ! strcasecmp(mySegid,segid) ) {
104 if ( resid == lastResid ) {
105 atomIndex[lastResid - firstResid + 1] = aid + 1;
106 }
else if ( resid < lastResid ) {
109 " out of order in segment " << segid <<
110 ", lookup for additional residues in this segment disabled.\n" <<
endi;
112 next->append(segid,resid,aid);
114 for ( ; lastResid < resid; ++lastResid ) atomIndex.add(aid);
115 atomIndex[lastResid - firstResid + 1] = aid + 1;
119 next->append(segid,resid,aid);
127 const char *segid,
int resid,
const char *aname)
const {
129 if (atomNames == NULL || resLookup == NULL)
131 NAMD_die(
"Tried to find atom from name on node other than node 0");
136 if ( resLookup->lookup(segid,resid,&i,&end) )
return -1;
137 for ( ; i < end; ++i ) {
138 #ifdef MEM_OPT_VERSION 139 Index idx = atomNames[i].atomnameIdx;
142 if ( ! strcasecmp(aname,atomNames[i].atomname) )
return i;
150 const char *segid,
int resid)
const {
152 if (atomNames == NULL || resLookup == NULL)
154 NAMD_die(
"Tried to find atom from name on node other than node 0");
158 if ( resLookup->lookup(segid,resid,&i,&end) )
return 0;
164 const char *segid,
int resid,
int index)
const {
166 if (atomNames == NULL || resLookup == NULL)
168 NAMD_die(
"Tried to find atom from name on node other than node 0");
172 if ( resLookup->lookup(segid,resid,&i,&end) )
return -1;
173 if ( index >= 0 && index < ( end - i ) )
return ( index + i );
190 if (
sizeof(
int32) != 4 ) {
NAMD_bug(
"sizeof(int32) != 4"); }
192 this->params = param;
200 is_lonepairs_psf = 1;
206 oneFourNbTholes=NULL;
211 lcpoParamType = NULL;
220 #ifdef MEM_OPT_VERSION 228 atomChargePool = NULL;
229 eachAtomCharge = NULL;
242 #ifndef MEM_OPT_VERSION 248 dihedralsByAtom=NULL;
249 impropersByAtom=NULL;
250 crosstermsByAtom=NULL;
252 gromacsPairByAtom=NULL;
257 oneFourNbTholesByAtom=NULL;
261 #ifdef MEM_OPT_VERSION 263 exclChkSigPool = NULL;
265 eachAtomExclSig = NULL;
267 fixedAtomsSet = NULL;
268 constrainedAtomsSet = NULL;
270 exclusionsByAtom=NULL;
271 fullExclusionsByAtom=NULL;
272 modExclusionsByAtom=NULL;
279 #ifdef MEM_OPT_VERSION 286 exPressureAtomFlags=NULL;
287 rigidBondLengths=NULL;
302 consTorqueIndexes=NULL;
303 consTorqueParams=NULL;
304 consForceIndexes=NULL;
307 moleculeStartIndex=NULL;
312 alch_unpert_bonds = NULL;
313 alch_unpert_angles = NULL;
314 alch_unpert_dihedrals = NULL;
324 #ifndef MEM_OPT_VERSION 352 numOneFourNbTholes=0;
363 numExPressureAtoms=0;
365 numFixedRigidBonds=0;
366 numMultipleDihedrals=0;
367 numMultipleImpropers=0;
376 numCalcFullExclusions=0;
377 numCalcOneFourNbTholes=0;
385 num_alch_unpert_Bonds = 0;
386 num_alch_unpert_Angles = 0;
387 num_alch_unpert_Dihedrals = 0;
461 #ifdef MEM_OPT_VERSION 463 read_mol_signatures(filename, param, cfgList);
465 read_psf_file(filename, param);
468 assignLCPOTypes( 0 );
482 #ifdef MEM_OPT_VERSION 483 NAMD_die(
"Sorry, plugin IO is not supported in the memory optimized version.");
487 int optflags = MOLFILE_BADOPTIONS;
488 molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*
sizeof(molfile_atom_t));
489 memset(atomarray, 0, natoms*
sizeof(molfile_atom_t));
492 int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
493 if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
495 NAMD_die(
"ERROR: plugin failed reading structure data");
497 if(optflags == MOLFILE_BADOPTIONS) {
499 NAMD_die(
"ERROR: plugin didn't initialize optional data flags");
501 if(optflags & MOLFILE_OCCUPANCY) {
502 setOccupancyData(atomarray);
504 if(optflags & MOLFILE_BFACTOR) {
505 setBFactorData(atomarray);
508 plgLoadAtomBasics(atomarray);
515 int *bondtype, nbondtypes;
517 if(pIOHdl->read_bonds!=NULL) {
518 if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder,
519 &bondtype, &nbondtypes, &bondtypename)){
520 NAMD_die(
"ERROR: failed reading bond information.");
525 plgLoadBonds(from,to);
529 int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
530 int ctermcols, ctermrows;
531 int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes;
532 int *impropertypes, numimpropertypes;
533 char **angletypenames, **dihedraltypenames, **impropertypenames;
535 plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
536 if(pIOHdl->read_angles!=NULL) {
537 if(pIOHdl->read_angles(pIOFileHdl,
538 &numAngles, &plgAngles,
539 &angletypes, &numangletypes, &angletypenames,
540 &numDihedrals, &plgDihedrals,
541 &dihedraltypes, &numdihedraltypes, &dihedraltypenames,
542 &numImpropers, &plgImpropers,
543 &impropertypes, &numimpropertypes, &impropertypenames,
544 &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) {
545 NAMD_die(
"ERROR: failed reading angle information.");
549 if(numAngles!=0) plgLoadAngles(plgAngles);
550 if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
551 if(numImpropers!=0) plgLoadImpropers(plgImpropers);
552 if(numCrossterms!=0) plgLoadCrossterms(plgCterms);
554 numRealBonds = numBonds;
558 assignLCPOTypes( 2 );
581 if (atomNames != NULL)
588 if (resLookup != NULL)
592 if (drudeConsts != NULL)
delete [] drudeConsts;
593 if (lphosts != NULL)
delete [] lphosts;
594 if (anisos != NULL)
delete [] anisos;
595 if (tholes != NULL)
delete [] tholes;
596 if (oneFourNbTholes != NULL)
delete[] oneFourNbTholes;
597 if (lphostIndexes != NULL)
delete [] lphostIndexes;
601 if (lcpoParamType != NULL)
delete [] lcpoParamType;
603 #ifdef MEM_OPT_VERSION 604 if(eachAtomSig)
delete [] eachAtomSig;
613 if (dihedrals != NULL)
616 if (impropers != NULL)
619 if (crossterms != NULL)
620 delete [] crossterms;
622 if (exclusions != NULL)
623 delete [] exclusions;
629 if (acceptors != NULL)
632 #ifdef MEM_OPT_VERSION 633 if(exclSigPool)
delete [] exclSigPool;
634 if(exclChkSigPool)
delete [] exclChkSigPool;
635 if(eachAtomExclSig)
delete [] eachAtomExclSig;
636 if(fixedAtomsSet)
delete fixedAtomsSet;
637 if(constrainedAtomsSet)
delete constrainedAtomsSet;
639 if (bondsByAtom != NULL)
640 delete [] bondsByAtom;
642 if (anglesByAtom != NULL)
643 delete [] anglesByAtom;
645 if (dihedralsByAtom != NULL)
646 delete [] dihedralsByAtom;
648 if (impropersByAtom != NULL)
649 delete [] impropersByAtom;
651 if (crosstermsByAtom != NULL)
652 delete [] crosstermsByAtom;
654 if (exclusionsByAtom != NULL)
655 delete [] exclusionsByAtom;
657 if (fullExclusionsByAtom != NULL)
658 delete [] fullExclusionsByAtom;
660 if (modExclusionsByAtom != NULL)
661 delete [] modExclusionsByAtom;
663 if (all_exclusions != NULL)
664 delete [] all_exclusions;
667 if (gromacsPairByAtom != NULL)
668 delete [] gromacsPairByAtom;
672 if (tholesByAtom != NULL)
673 delete [] tholesByAtom;
674 if (anisosByAtom != NULL)
675 delete [] anisosByAtom;
676 if (oneFourNbTholesByAtom != NULL)
677 delete[] oneFourNbTholesByAtom;
682 if (lcpoParamType != NULL)
683 delete [] lcpoParamType;
685 if (fixedAtomFlags != NULL)
686 delete [] fixedAtomFlags;
688 if (stirIndexes != NULL)
689 delete [] stirIndexes;
692 #ifdef MEM_OPT_VERSION 693 if(clusterSigs != NULL){
694 delete [] clusterSigs;
700 if (clusterSize != NULL)
701 delete [] clusterSize;
703 if (exPressureAtomFlags != NULL)
704 delete [] exPressureAtomFlags;
706 if (rigidBondLengths != NULL)
707 delete [] rigidBondLengths;
710 if (fepAtomFlags != NULL)
711 delete [] fepAtomFlags;
712 if (alch_unpert_bonds != NULL)
713 delete [] alch_unpert_bonds;
714 if (alch_unpert_angles != NULL)
715 delete [] alch_unpert_angles;
716 if (alch_unpert_dihedrals != NULL)
717 delete [] alch_unpert_dihedrals;
721 if (ssAtomFlags != NULL)
722 delete [] ssAtomFlags;
723 if (ss_vdw_type != NULL)
724 delete [] ss_vdw_type;
725 if (ss_index != NULL)
729 if (qmAtomGroup != NULL)
730 delete [] qmAtomGroup;
732 if (qmAtmIndx != NULL)
735 if (qmAtmChrg != NULL)
739 if (qmGrpNumBonds != NULL)
740 delete [] qmGrpNumBonds;
742 if (qmGrpSizes != NULL)
743 delete [] qmGrpSizes;
745 if (qmDummyBondVal != NULL)
746 delete [] qmDummyBondVal;
748 if (qmMMNumTargs != NULL)
749 delete [] qmMMNumTargs;
754 if (qmGrpChrg != NULL)
757 if (qmGrpMult != NULL)
760 if (qmMeMMindx != NULL)
761 delete [] qmMeMMindx;
763 if (qmMeQMGrp != NULL)
766 if (qmLSSSize != NULL)
769 if (qmLSSIdxs != NULL)
772 if (qmLSSMass != NULL)
775 if (qmLSSRefSize != NULL)
776 delete [] qmLSSRefSize;
778 if (qmLSSRefIDs != NULL)
779 delete [] qmLSSRefIDs;
781 if (qmLSSRefMass != NULL)
782 delete [] qmLSSRefMass;
784 if (qmMMBond != NULL) {
785 for(
int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
786 if (qmMMBond[grpIndx] != NULL)
787 delete [] qmMMBond[grpIndx];
792 if (qmGrpBonds != NULL) {
793 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
794 if (qmGrpBonds[grpIndx] != NULL)
795 delete [] qmGrpBonds[grpIndx];
797 delete [] qmGrpBonds;
800 if (qmMMBondedIndx != NULL) {
801 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
802 if (qmMMBondedIndx[grpIndx] != NULL)
803 delete [] qmMMBondedIndx[grpIndx];
805 delete [] qmMMBondedIndx;
808 if (qmMMChargeTarget != NULL) {
809 for(
int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
810 if (qmMMChargeTarget[grpIndx] != NULL)
811 delete [] qmMMChargeTarget[grpIndx];
813 delete [] qmMMChargeTarget;
816 if (qmElementArray != NULL){
817 for(
int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
818 if (qmElementArray[atmInd] != NULL)
819 delete [] qmElementArray[atmInd];
821 delete [] qmElementArray;
824 if (qmDummyElement != NULL){
825 for(
int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
826 if (qmDummyElement[atmInd] != NULL)
827 delete [] qmDummyElement[atmInd];
829 delete [] qmDummyElement;
832 if (qmCustPCSizes != NULL){
833 delete [] qmCustPCSizes;
836 if (qmCustomPCIdxs != NULL){
837 delete [] qmCustomPCIdxs;
840 if (cSMDindex != NULL) {
841 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
842 if (cSMDindex[grpIndx] != NULL)
843 delete [] cSMDindex[grpIndx];
848 if (cSMDpairs != NULL) {
849 for(
int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) {
850 if (cSMDpairs[instIndx] != NULL)
851 delete [] cSMDpairs[instIndx];
856 if (cSMDindxLen != NULL)
857 delete [] cSMDindxLen;
860 if (cSMDVels != NULL)
862 if (cSMDcoffs != NULL)
865 #ifndef MEM_OPT_VERSION 872 #ifndef MEM_OPT_VERSION 893 void Molecule::read_psf_file(
char *fname,
Parameters *params)
903 if ( (psf_file =
Fopen(fname,
"r")) == NULL)
905 sprintf(err_msg,
"UNABLE TO OPEN .psf FILE %s", fname);
921 sprintf(err_msg,
"EMPTY .psf FILE %s", fname);
929 sprintf(err_msg,
"UNABLE TO FIND \"PSF\" STRING IN PSF FILE %s",
938 iout <<
iWARN <<
"Reading PSF supporting DRUDE without " 939 "enabling the Drude model in the simulation config file\n" <<
endi;
958 sprintf(err_msg,
"MISSING EVERYTHING BUT PSF FROM %s", fname);
966 sprintf(err_msg,
"CAN NOT FIND \"NTITLE\" STRING IN PSF FILE %s",
971 sscanf(buffer,
"%d", &NumTitle);
987 sprintf(err_msg,
"FOUND EOF INSTEAD OF NATOM IN PSF FILE %s",
1000 sprintf(err_msg,
"DIDN'T FIND \"NATOM\" IN PSF FILE %s",
1006 sscanf(buffer,
"%d", &numAtoms);
1008 read_atoms(psf_file, params);
1021 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NBONDS IN PSF");
1027 NAMD_die(
"DID NOT FIND NBOND AFTER ATOM LIST IN PSF");
1031 sscanf(buffer,
"%d", &numBonds);
1034 read_bonds(psf_file);
1047 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NTHETA IN PSF");
1053 NAMD_die(
"DID NOT FIND NTHETA AFTER BOND LIST IN PSF");
1057 sscanf(buffer,
"%d", &numAngles);
1060 read_angles(psf_file, params);
1073 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NPHI IN PSF");
1079 NAMD_die(
"DID NOT FIND NPHI AFTER ANGLE LIST IN PSF");
1083 sscanf(buffer,
"%d", &numDihedrals);
1086 read_dihedrals(psf_file, params);
1099 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NIMPHI IN PSF");
1105 NAMD_die(
"DID NOT FIND NIMPHI AFTER ATOM LIST IN PSF");
1109 sscanf(buffer,
"%d", &numImpropers);
1112 read_impropers(psf_file, params);
1125 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NDON IN PSF");
1131 NAMD_die(
"DID NOT FIND NDON AFTER ATOM LIST IN PSF");
1135 sscanf(buffer,
"%d", &numDonors);
1138 read_donors(psf_file);
1151 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NACC IN PSF");
1157 NAMD_die(
"DID NOT FIND NACC AFTER ATOM LIST IN PSF");
1161 sscanf(buffer,
"%d", &numAcceptors);
1164 read_acceptors(psf_file);
1173 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NNB IN PSF FILE");
1178 sscanf(buffer,
"%d", &numExclusions);
1181 read_exclusions(psf_file);
1200 int is_found_numlp = 0;
1201 int is_found_numaniso = 0;
1202 int is_found_ncrterm = 0;
1203 while ( ! is_found ) {
1209 if (ret_code != 0) {
1212 else if ( (is_found_numlp =
NAMD_find_word(buffer,
"NUMLP")) > 0) {
1213 is_found = is_found_numlp;
1215 else if ( (is_found_numaniso =
NAMD_find_word(buffer,
"NUMANISO")) > 0) {
1216 is_found = is_found_numaniso;
1218 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1219 is_found = is_found_ncrterm;
1223 if (is_found_numlp) {
1226 sscanf(buffer,
"%d", &numLphosts);
1229 if (numLphosts == 0) {
1234 is_lonepairs_psf = 0;
1235 }
else if (
simParams->CUDASOAintegrate) {
1243 NAMD_die(
"FOUND LONE PAIR HOSTS IN PSF WITH \"LONEPAIRS\" DISABLED IN CONFIG FILE");
1247 if (numBonds) process_bonds(params);
1249 if (is_found_numlp) {
1251 if (numLphosts > 0) read_lphosts(psf_file);
1255 is_found_numaniso = 0;
1256 is_found_ncrterm = 0;
1257 while ( ! is_found ) {
1263 if (ret_code != 0) {
1266 else if ( (is_found_numaniso =
NAMD_find_word(buffer,
"NUMANISO")) > 0) {
1267 is_found = is_found_numaniso;
1269 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1270 is_found = is_found_ncrterm;
1275 if (is_found_numaniso && is_drude_psf) {
1278 sscanf(buffer,
"%d", &numAnisos);
1279 if (numAnisos > 0) read_anisos(psf_file);
1283 is_found_ncrterm = 0;
1284 while ( ! is_found ) {
1290 if (ret_code != 0) {
1293 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1294 is_found = is_found_ncrterm;
1298 else if (is_drude_psf ) {
1299 NAMD_die(
"DID NOT FIND REQUIRED NUMANISO IN DRUDE PSF FILE");
1301 else if (is_found_numaniso ) {
1302 NAMD_die(
"FOUND NUMANISO IN PSF FILE MISSING DRUDE DESIGNATION");
1305 if (is_found_ncrterm) {
1308 sscanf(buffer,
"%d", &numCrossterms);
1309 if (numCrossterms > 0) read_crossterms(psf_file, params);
1317 if (is_drude_psf && numDrudeAtoms) {
1319 Bond *newbonds =
new Bond[numBonds+numDrudeAtoms];
1320 memcpy(newbonds, bonds, numBonds*
sizeof(
Bond));
1323 int origNumBonds = numBonds;
1324 for (i=0; i < numAtoms; i++) {
1325 if (!is_drude(i))
continue;
1326 Bond *b = &(bonds[numBonds++]);
1330 atomNames[i-1].atomtype, atomNames[i].atomtype, b
1333 if (numBonds-origNumBonds != numDrudeAtoms) {
1334 NAMD_die(
"must have same number of Drude particles and parents");
1339 numRealBonds = numBonds;
1340 build_atom_status();
1361 void Molecule::read_atoms(FILE *fd,
Parameters *params)
1366 int last_atom_number=0;
1368 char segment_name[11];
1369 char residue_number[11];
1370 char residue_name[11];
1378 atoms =
new Atom[numAtoms];
1390 hydrogenGroup.resize(0);
1392 if (atoms == NULL || atomNames == NULL )
1394 NAMD_die(
"memory allocation failed in Molecule::read_atoms");
1400 while (atom_number < numAtoms)
1413 read_count=sscanf(buffer,
"%d %s %s %s %s %s %f %f",
1414 &atom_number, segment_name, residue_number,
1415 residue_name, atom_name, atom_type, &charge, &mass);
1416 if (mass <= 0.05) ++numZeroMassAtoms;
1419 if (read_count != 8)
1423 sprintf(err_msg,
"BAD ATOM LINE FORMAT IN PSF FILE IN ATOM LINE %d\nLINE=%s",
1424 last_atom_number+1, buffer);
1437 read_count=sscanf(buffer,
1440 "%*d %*s %*s %*s %*s %*s %*f %*f %*d %f %f", &alpha, &
thole);
1441 if (read_count != 2)
1445 sprintf(err_msg,
"BAD ATOM LINE FORMAT IN PSF FILE " 1446 "IN ATOM LINE %d\nLINE=%s", last_atom_number+1, buffer);
1449 drudeConsts[atom_number-1].alpha = alpha;
1450 drudeConsts[atom_number-1].thole =
thole;
1451 if (fabs(alpha) >= 1e-6) ++numDrudeAtoms;
1457 if ( sscanf(atom_type,
"%d", &atom_type_num) > 0 )
1459 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.");
1463 if (atom_number != last_atom_number+1)
1467 sprintf(err_msg,
"ATOM NUMBERS OUT OF ORDER AT ATOM #%d OF PSF FILE",
1468 last_atom_number+1);
1477 int reslength = strlen(residue_name)+1;
1478 int namelength = strlen(atom_name)+1;
1479 int typelength = strlen(atom_type)+1;
1481 atomNames[atom_number-1].resname = nameArena->getNewArray(reslength);
1482 atomNames[atom_number-1].atomname = nameArena->getNewArray(namelength);
1483 atomNames[atom_number-1].atomtype = nameArena->getNewArray(typelength);
1485 if (atomNames[atom_number-1].resname == NULL)
1487 NAMD_die(
"memory allocation failed in Molecule::read_atoms");
1491 strcpy(atomNames[atom_number-1].resname, residue_name);
1492 strcpy(atomNames[atom_number-1].atomname, atom_name);
1493 strcpy(atomNames[atom_number-1].atomtype, atom_type);
1494 atoms[atom_number-1].mass = mass;
1495 atoms[atom_number-1].charge = charge;
1499 if ( tmpResLookup ) tmpResLookup =
1500 tmpResLookup->
append(segment_name, atoi(residue_number), atom_number-1);
1504 memcpy(one->
segname, segment_name, strlen(segment_name)+1);
1505 one->
resid = atoi(residue_number);
1510 }
else if (atoms[atom_number-1].mass <= 0.05) {
1512 }
else if (atoms[atom_number-1].mass < 1.0) {
1513 atoms[atom_number-1].status |=
DrudeAtom;
1514 }
else if (atoms[atom_number-1].mass <=3.5) {
1516 }
else if ((atomNames[atom_number-1].atomname[0] ==
'O') &&
1517 (atoms[atom_number-1].mass >= 14.0) &&
1518 (atoms[atom_number-1].mass <= 18.0)) {
1524 &(atoms[atom_number-1]));
1543 void Molecule::read_bonds(FILE *fd)
1551 bonds=
new Bond[numBonds];
1555 NAMD_die(
"memory allocations failed in Molecule::read_bonds");
1559 while (num_read < numBonds)
1571 if (atom_nums[j] >= numAtoms)
1575 sprintf(err_msg,
"BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1581 Bond *b = &(bonds[num_read]);
1582 b->
atom1=atom_nums[0];
1583 b->
atom2=atom_nums[1];
1592 void Molecule::process_bonds(
Parameters *params) {
1595 int origNumBonds = numBonds;
1597 int numZeroFrcBonds = 0;
1599 int numDrudeBonds = 0;
1605 bool bond_found =
false;
1606 for (
int old_read=0; old_read < origNumBonds; ++old_read)
1609 Bond *b = &(bonds[num_read]);
1610 *b = bonds[old_read];
1611 atom_nums[0] = b->
atom1;
1612 atom_nums[1] = b->
atom2;
1617 atomNames[atom_nums[0]].atomtype,
1618 atomNames[atom_nums[1]].atomtype,
1632 if (!is_lp_bond && !bond_found) {
1634 snprintf(err_msg,
sizeof(err_msg),
1635 "UNABLE TO FIND BOND PARAMETERS FOR %s %s (ATOMS %i %i)",
1636 atomNames[atom_nums[0]].atomtype, atomNames[atom_nums[1]].atomtype,
1641 numZeroFrcBonds += (k == 0.);
1642 numLPBonds += is_lp_bond;
1643 numDrudeBonds += is_drude_bond;
1653 if ( is_lonepairs_psf ) {
1654 if (k == 0. || is_lp_bond) --numBonds;
1657 numLPBonds -= is_lp_bond;
1658 if (k == 0. && !is_lp_bond) --numBonds;
1667 if (k == 0. || is_lp_bond || is_drude_bond) --numBonds;
1673 if ( num_read != numBonds ) {
1674 NAMD_bug(
"num_read != numBonds in Molecule::process_bonds()");
1678 if ( numBonds != origNumBonds ) {
1679 if (numZeroFrcBonds) {
1680 iout <<
iWARN <<
"Ignored " << numZeroFrcBonds <<
1681 " bonds with zero force constants.\n" <<
1682 iWARN <<
"Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
1686 iout <<
iWARN <<
"Ignored " << numLPBonds <<
1687 " bonds with lone pairs.\n" <<
1688 iWARN <<
"Will infer lonepair bonds from LPhost entries.\n" <<
endi;
1690 if (numDrudeBonds) {
1691 iout <<
iWARN <<
"Ignored " << numDrudeBonds <<
1692 " bonds with Drude particles.\n" <<
1693 iWARN <<
"Will use polarizability to assign Drude bonds.\n" <<
endi;
1706 alch_unpert_bonds=
new Bond[num_alch_unpert_Bonds];
1708 if (alch_unpert_bonds == NULL) {
1709 NAMD_die(
"memory allocations failed in Molecule::read_alch_unpert_bonds");
1712 while (num_read < num_alch_unpert_Bonds) {
1713 for (j=0; j<2; j++) {
1716 if (atom_nums[j] >= numAtoms) {
1719 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);
1724 Bond *b = &(alch_unpert_bonds[num_read]);
1725 b->
atom1=atom_nums[0];
1726 b->
atom2=atom_nums[1];
1750 void Molecule::read_angles(FILE *fd,
Parameters *params)
1756 int origNumAngles = numAngles;
1758 angles=
new Angle[numAngles];
1762 NAMD_die(
"memory allocation failed in Molecule::read_angles");
1766 while (num_read < numAngles)
1779 if (atom_nums[j] >= numAtoms)
1783 sprintf(err_msg,
"ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1789 angles[num_read].atom1=atom_nums[0];
1790 angles[num_read].atom2=atom_nums[1];
1791 angles[num_read].atom3=atom_nums[2];
1796 atomNames[atom_nums[0]].atomtype,
1797 atomNames[atom_nums[1]].atomtype,
1798 atomNames[atom_nums[2]].atomtype,
1799 &(angles[num_read]),
simParams->alchOn ? -1 : 0);
1800 if ( angles[num_read].angle_type == -1 ) {
1801 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n" 1806 Real k, t0, k_ub, r_ub;
1807 if ( angles[num_read].angle_type == -1 ) { k = -1.; k_ub = -1.; }
else 1809 if ( k == 0. && k_ub == 0. ) --numAngles;
1814 if ( numAngles != origNumAngles ) {
1815 iout <<
iWARN <<
"Ignored " << origNumAngles - numAngles <<
1816 " angles with zero force constants.\n" <<
endi;
1829 alch_unpert_angles=
new Angle[num_alch_unpert_Angles];
1831 if (alch_unpert_angles == NULL) {
1832 NAMD_die(
"memory allocation failed in Molecule::read_alch_unpert_angles");
1835 while (num_read < num_alch_unpert_Angles) {
1836 for (j=0; j<3; j++) {
1839 if (atom_nums[j] >= numAtoms) {
1841 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);
1846 alch_unpert_angles[num_read].atom1=atom_nums[0];
1847 alch_unpert_angles[num_read].atom2=atom_nums[1];
1848 alch_unpert_angles[num_read].atom3=atom_nums[2];
1871 void Molecule::read_dihedrals(FILE *fd,
Parameters *params)
1874 int last_atom_nums[4];
1878 Bool duplicate_bond;
1883 last_atom_nums[j] = -1;
1886 dihedrals =
new Dihedral[numDihedrals];
1888 if (dihedrals == NULL)
1890 NAMD_die(
"memory allocation failed in Molecule::read_dihedrals");
1894 while (num_read < numDihedrals)
1896 duplicate_bond =
TRUE;
1908 if (atom_nums[j] >= numAtoms)
1912 sprintf(err_msg,
"DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1917 if (atom_nums[j] != last_atom_nums[j])
1919 duplicate_bond =
FALSE;
1922 last_atom_nums[j] = atom_nums[j];
1932 if (multiplicity == 2)
1934 numMultipleDihedrals++;
1944 dihedrals[num_unique-1].atom1=atom_nums[0];
1945 dihedrals[num_unique-1].atom2=atom_nums[1];
1946 dihedrals[num_unique-1].atom3=atom_nums[2];
1947 dihedrals[num_unique-1].atom4=atom_nums[3];
1951 atomNames[atom_nums[0]].atomtype,
1952 atomNames[atom_nums[1]].atomtype,
1953 atomNames[atom_nums[2]].atomtype,
1954 atomNames[atom_nums[3]].atomtype,
1955 &(dihedrals[num_unique-1]),
1956 multiplicity,
simParams->alchOn ? -1 : 0);
1957 if ( dihedrals[num_unique-1].dihedral_type == -1 ) {
1958 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n" 1965 numDihedrals = num_unique;
1977 alch_unpert_dihedrals =
new Dihedral[num_alch_unpert_Dihedrals];
1979 if (alch_unpert_dihedrals == NULL) {
1980 NAMD_die(
"memory allocation failed in Molecule::read_alch_unpert_dihedrals");
1983 while (num_read < num_alch_unpert_Dihedrals) {
1984 for (j=0; j<4; j++) {
1987 if (atom_nums[j] >= numAtoms) {
1990 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);
1995 alch_unpert_dihedrals[num_read].atom1=atom_nums[0];
1996 alch_unpert_dihedrals[num_read].atom2=atom_nums[1];
1997 alch_unpert_dihedrals[num_read].atom3=atom_nums[2];
1998 alch_unpert_dihedrals[num_read].atom4=atom_nums[3];
2021 void Molecule::read_impropers(FILE *fd,
Parameters *params)
2024 int last_atom_nums[4];
2028 Bool duplicate_bond;
2034 last_atom_nums[j] = -1;
2037 impropers=
new Improper[numImpropers];
2039 if (impropers == NULL)
2041 NAMD_die(
"memory allocation failed in Molecule::read_impropers");
2045 while (num_read < numImpropers)
2047 duplicate_bond =
TRUE;
2059 if (atom_nums[j] >= numAtoms)
2063 sprintf(err_msg,
"IMPROPERS INDEX %d GREATER THAN NATOM %d IN IMPROPERS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
2067 if (atom_nums[j] != last_atom_nums[j])
2069 duplicate_bond =
FALSE;
2072 last_atom_nums[j] = atom_nums[j];
2083 if (multiplicity == 2)
2086 numMultipleImpropers++;
2097 impropers[num_unique-1].atom1=atom_nums[0];
2098 impropers[num_unique-1].atom2=atom_nums[1];
2099 impropers[num_unique-1].atom3=atom_nums[2];
2100 impropers[num_unique-1].atom4=atom_nums[3];
2104 atomNames[atom_nums[0]].atomtype,
2105 atomNames[atom_nums[1]].atomtype,
2106 atomNames[atom_nums[2]].atomtype,
2107 atomNames[atom_nums[3]].atomtype,
2108 &(impropers[num_unique-1]),
2117 numImpropers = num_unique;
2137 void Molecule::read_crossterms(FILE *fd,
Parameters *params)
2141 int last_atom_nums[8];
2144 Bool duplicate_bond;
2149 last_atom_nums[j] = -1;
2152 crossterms=
new Crossterm[numCrossterms];
2154 if (crossterms == NULL)
2156 NAMD_die(
"memory allocation failed in Molecule::read_crossterms");
2160 while (num_read < numCrossterms)
2162 duplicate_bond =
TRUE;
2174 if (atom_nums[j] >= numAtoms)
2178 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);
2182 if (atom_nums[j] != last_atom_nums[j])
2184 duplicate_bond =
FALSE;
2187 last_atom_nums[j] = atom_nums[j];
2193 iout <<
iWARN <<
"Duplicate cross-term detected.\n" <<
endi;
2197 crossterms[num_read].atom1=atom_nums[0];
2198 crossterms[num_read].atom2=atom_nums[1];
2199 crossterms[num_read].atom3=atom_nums[2];
2200 crossterms[num_read].atom4=atom_nums[3];
2201 crossterms[num_read].atom5=atom_nums[4];
2202 crossterms[num_read].atom6=atom_nums[5];
2203 crossterms[num_read].atom7=atom_nums[6];
2204 crossterms[num_read].atom8=atom_nums[7];
2208 atomNames[atom_nums[0]].atomtype,
2209 atomNames[atom_nums[1]].atomtype,
2210 atomNames[atom_nums[2]].atomtype,
2211 atomNames[atom_nums[3]].atomtype,
2212 atomNames[atom_nums[4]].atomtype,
2213 atomNames[atom_nums[5]].atomtype,
2214 atomNames[atom_nums[6]].atomtype,
2215 atomNames[atom_nums[7]].atomtype,
2216 &(crossterms[num_read]));
2218 if(!duplicate_bond) num_read++;
2221 numCrossterms = num_read;
2244 void Molecule::read_donors(FILE *fd)
2252 donors=
new Bond[numDonors];
2256 NAMD_die(
"memory allocations failed in Molecule::read_donors");
2260 while (num_read < numDonors)
2272 if (d[j] >= numAtoms)
2277 "DONOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE",
2278 d[j]+1, numAtoms, num_read+1);
2288 Bond *b = &(donors[num_read]);
2318 void Molecule::read_acceptors(FILE *fd)
2326 acceptors=
new Bond[numAcceptors];
2328 if (acceptors == NULL)
2330 NAMD_die(
"memory allocations failed in Molecule::read_acceptors");
2334 while (num_read < numAcceptors)
2346 if (d[j] >= numAtoms)
2350 sprintf(err_msg,
"ACCEPTOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE", d[j]+1, numAtoms, num_read+1);
2360 Bond *b = &(acceptors[num_read]);
2399 void Molecule::read_exclusions(FILE *fd)
2402 int *exclusion_atoms;
2403 register int num_read=0;
2406 register int insert_index=0;
2410 exclusions =
new Exclusion[numExclusions];
2411 exclusion_atoms =
new int[numExclusions];
2413 if ( (exclusions == NULL) || (exclusion_atoms == NULL) )
2415 NAMD_die(
"memory allocation failed in Molecule::read_exclusions");
2419 for (num_read=0; num_read<numExclusions; num_read++)
2427 if (exclusion_atoms[num_read] >= numAtoms)
2431 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);
2440 for (num_read=0; num_read<numAtoms; num_read++)
2446 if (current_index>numExclusions)
2450 sprintf(err_msg,
"EXCLUSION INDEX %d LARGER THAN NUMBER OF EXLCUSIONS %d IN PSF FILE, EXCLUSION #%d\n",
2451 current_index+1, numExclusions, num_read);
2458 if (current_index != last_index)
2464 for (insert_index=last_index;
2465 insert_index<current_index; insert_index++)
2472 int a2 = exclusion_atoms[insert_index];
2474 exclusions[insert_index].atom1 = a1;
2475 exclusions[insert_index].atom2 = a2;
2476 }
else if ( a2 < a1 ) {
2477 exclusions[insert_index].atom1 = a2;
2478 exclusions[insert_index].atom2 = a1;
2481 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN PSF FILE\n", a1+1);
2486 last_index=current_index;
2491 delete [] exclusion_atoms;
2507 void Molecule::read_exclusions(
int* atom_i,
int* atom_j,
int num_exclusion)
2511 exclusions =
new Exclusion[num_exclusion];
2512 int loop_counter = 0;
2516 if ( (exclusions == NULL) )
2518 NAMD_die(
"memory allocation failed in Molecule::read_exclusions");
2522 for (loop_counter = 0; loop_counter < num_exclusion; loop_counter++) {
2524 if ( (atom_i == NULL) || (atom_j == NULL) ) {
2525 NAMD_die(
"null pointer expection in Molecule::read_exclusions");
2528 a = atom_i[loop_counter];
2529 b = atom_j[loop_counter];
2531 exclusions[loop_counter].atom1 = a;
2532 exclusions[loop_counter].atom2 = b;
2534 exclusions[loop_counter].atom1 = b;
2535 exclusions[loop_counter].atom2 = a;
2537 exclusionSet.add(
Exclusion(exclusions[loop_counter].atom1,exclusions[loop_counter].atom2));
2541 iout <<
iINFO <<
"ADDED " << num_exclusion <<
" EXPLICIT EXCLUSIONS: THIS VALUE WILL *NOT* BE ADDED TO THE STRUCTURE SUMMARY\n" <<
endi;
2565 void Molecule::read_lphosts(FILE *fd)
2577 Real value1, value2, value3;
2580 lphosts =
new Lphost[numLphosts];
2581 if (lphosts == NULL) {
2582 NAMD_die(
"memory allocation failed in Molecule::read_lphosts");
2584 for (i = 0; i < numLphosts; i++) {
2587 read_count=sscanf(buffer,
"%d %d %6s %f %f %f",
2588 &numhosts, &index, weight, &value1, &value2, &value3);
2593 if (read_count != 6 || index != next_index ||
2594 strcmp(weight,
"F") != 0 || numhosts < 2 || 3 < numhosts) {
2596 sprintf(err_msg,
"BAD FORMAT FOR LPHOST LINE %d IN PSF FILE LINE\n" 2597 "LINE=%s\n", i+1, buffer);
2600 lphosts[i].numhosts = numhosts;
2601 next_index += numhosts + 1;
2602 if (numhosts == 2) {
2603 lphosts[i].distance = value1;
2604 lphosts[i].angle = value2;
2605 lphosts[i].dihedral = 0.0;
2608 lphosts[i].distance = value1;
2609 lphosts[i].angle = value2 * (
M_PI/180);
2610 lphosts[i].dihedral = value3 * (
M_PI/180);
2614 Bond *newbonds =
new Bond[numBonds+numLphosts];
2615 memcpy(newbonds, bonds, numBonds*
sizeof(
Bond));
2619 params->add_bond_param(
"LP X 0.0 0.0\n",
FALSE);
2620 const char dummyLPtype[] =
"LP";
2621 const char dummyHostType[] =
"X";
2622 for (i = 0; i < numLphosts; i++) {
2628 lphosts[i].atom4 = ( lphosts[i].numhosts == 3 ?
2631 Bond *b = &(bonds[numBonds++]);
2632 b->
atom1 = lphosts[i].atom1;
2633 b->
atom2 = lphosts[i].atom2;
2649 void Molecule::read_anisos(FILE *fd)
2652 int numhosts, index, i, read_count;
2655 anisos =
new Aniso[numAnisos];
2658 NAMD_die(
"memory allocation failed in Molecule::read_anisos");
2660 for (i = 0; i < numAnisos; i++)
2664 read_count=sscanf(buffer,
"%f %f %f", &k11, &k22, &k33);
2665 if (read_count != 3)
2668 sprintf(err_msg,
"BAD FORMAT FOR ANISO LINE %d IN PSF FILE LINE\n" 2669 "LINE=%s\n", i+1, buffer);
2672 anisos[i].k11 = k11;
2673 anisos[i].k22 = k22;
2674 anisos[i].k33 = k33;
2676 for (i = 0; i < numAnisos; i++) {
2689 if (atomType[0] ==
'H' || atomType[0] ==
'h') {
2693 }
else if (atomType[0] ==
'C' || atomType[0] ==
'c') {
2696 strcmp(atomType,
"CT" )==0 )
2704 else if (numBonds == 2)
2706 else if (numBonds == 3)
2708 else if (numBonds == 4)
2716 else if (numBonds == 3)
2723 }
else if (atomType[0] ==
'N' || atomType[0] ==
'n') {
2724 if ( strcmp(atomType,
"N3" ) == 0 ) {
2727 else if (numBonds == 2)
2729 else if (numBonds == 3)
2737 else if (numBonds == 2)
2739 else if (numBonds == 3)
2746 }
else if (atomType[0] ==
'O' || atomType[0] ==
'o') {
2748 if ( strcmp(atomType,
"O" )==0) {
2750 }
else if (strcmp(atomType,
"O2" )==0) {
2755 else if (numBonds == 2)
2762 }
else if (atomType[0] ==
'S' || atomType[0] ==
's') {
2763 if ( strcmp(atomType,
"SH" )==0) {
2770 }
else if (atomType[0] ==
'P' || atomType[0] ==
'p') {
2773 else if (numBonds == 4)
2777 }
else if (atomType[0] ==
'Z') {
2779 }
else if ( strcmp(atomType,
"MG" )==0) {
2790 if (atomType[0] ==
'H') {
2794 }
else if (atomType[0] ==
'C') {
2796 atomType[1] ==
'T' ||
2797 strcmp(atomType,
"CP1" )==0 ||
2798 strcmp(atomType,
"CP2" )==0 ||
2799 strcmp(atomType,
"CP3" )==0 ||
2800 strcmp(atomType,
"CS" )==0 ) {
2803 else if (numBonds == 2)
2805 else if (numBonds == 3)
2807 else if (numBonds == 4)
2813 strcmp(atomType,
"C" )==0 ||
2814 strcmp(atomType,
"CA" )==0 ||
2815 strcmp(atomType,
"CC" )==0 ||
2816 strcmp(atomType,
"CD" )==0 ||
2817 strcmp(atomType,
"CN" )==0 ||
2818 strcmp(atomType,
"CY" )==0 ||
2819 strcmp(atomType,
"C3" )==0 ||
2820 strcmp(atomType,
"CE1" )==0 ||
2821 strcmp(atomType,
"CE2" )==0 ||
2822 strcmp(atomType,
"CST" )==0 ||
2823 strcmp(atomType,
"CAP" )==0 ||
2824 strcmp(atomType,
"COA" )==0 ||
2825 strcmp(atomType,
"CPT" )==0 ||
2826 strcmp(atomType,
"CPH1")==0 ||
2827 strcmp(atomType,
"CPH2")==0
2831 else if (numBonds == 3)
2840 }
else if (atomType[0] ==
'N') {
2845 strcmp(atomType,
"NH3" )==0 ||
2848 strcmp(atomType,
"NP" )==0
2852 else if (numBonds == 2)
2854 else if (numBonds == 3)
2860 strcmp(atomType,
"NY" )==0 ||
2861 strcmp(atomType,
"NC2" )==0 ||
2862 strcmp(atomType,
"N" )==0 ||
2863 strcmp(atomType,
"NH1" )==0 ||
2864 strcmp(atomType,
"NH2" )==0 ||
2865 strcmp(atomType,
"NR1" )==0 ||
2866 strcmp(atomType,
"NR2" )==0 ||
2867 strcmp(atomType,
"NR3" )==0 ||
2868 strcmp(atomType,
"NPH" )==0 ||
2869 strcmp(atomType,
"NC" )==0
2873 else if (numBonds == 2)
2875 else if (numBonds == 3)
2884 }
else if (atomType[0] ==
'O') {
2886 strcmp(atomType,
"OH1" )==0 ||
2887 strcmp(atomType,
"OS" )==0 ||
2888 strcmp(atomType,
"OC" )==0 ||
2889 strcmp(atomType,
"OT" )==0
2893 else if (numBonds == 2)
2898 strcmp(atomType,
"O" )==0 ||
2899 strcmp(atomType,
"OB" )==0 ||
2900 strcmp(atomType,
"OST" )==0 ||
2901 strcmp(atomType,
"OCA" )==0 ||
2902 strcmp(atomType,
"OM" )==0
2906 strcmp(atomType,
"OC" )==0
2914 }
else if (atomType[0] ==
'S') {
2921 }
else if (atomType[0] ==
'P') {
2924 else if (numBonds == 4)
2939 void Molecule::assignLCPOTypes(
int inputType) {
2940 int *heavyBonds =
new int[numAtoms];
2941 for (
int i = 0; i < numAtoms; i++)
2943 for (
int i = 0; i < numBonds; i++ ) {
2944 Bond curBond = bonds[i];
2945 int a1 = bonds[i].
atom1;
2946 int a2 = bonds[i].atom2;
2947 if (atoms[a1].mass > 2.f && atoms[a2].mass > 2.f) {
2953 lcpoParamType =
new int[numAtoms];
2956 for (
int i = 0; i < numAtoms; i++) {
2959 if (inputType == 1) {
2973 if ( atoms[i].mass < 1.5 && lcpoParamType[i] != 0 ) {
2976 lcpoParamType[i] = 0;
2980 lcpoParamType[i] = 0;
2983 CkPrintf(
"ERROR in Molecule::assignLCPOTypes(): " 2984 "Light atom given heavy atom LCPO type.\n");
2992 iout <<
iWARN <<
"LONE PAIRS TO BE IGNORED BY SASA\n" <<
endi;
2995 iout <<
iWARN <<
"DRUDE PARTICLES TO BE IGNORED BY SASA\n" <<
endi;
2998 delete [] heavyBonds;
3002 void Molecule::plgLoadAtomBasics(molfile_atom_t *atomarray){
3003 atoms =
new Atom[numAtoms];
3008 hydrogenGroup.resize(0);
3012 for(
int i=0; i<numAtoms; i++) {
3013 int reslength = strlen(atomarray[i].resname)+1;
3014 int namelength = strlen(atomarray[i].name)+1;
3015 int typelength = strlen(atomarray[i].type)+1;
3016 atomNames[i].resname = nameArena->getNewArray(reslength);
3017 atomNames[i].atomname = nameArena->getNewArray(namelength);
3018 atomNames[i].atomtype = nameArena->getNewArray(typelength);
3019 strcpy(atomNames[i].resname, atomarray[i].resname);
3020 strcpy(atomNames[i].atomname, atomarray[i].name);
3021 strcpy(atomNames[i].atomtype, atomarray[i].type);
3023 atoms[i].mass = atomarray[i].mass;
3024 atoms[i].charge = atomarray[i].charge;
3029 tmpResLookup = tmpResLookup->
append(atomarray[i].segid, atomarray[i].resid, i);
3034 memcpy(one->
segname, atomarray[i].segid, strlen(atomarray[i].segid)+1);
3035 one->
resid = atomarray[i].resid;
3039 }
else if(atoms[i].mass <= 0.05) {
3041 }
else if(atoms[i].mass < 1.0) {
3043 }
else if(atoms[i].mass <= 3.5) {
3045 }
else if((atomNames[i].atomname[0] ==
'O') &&
3046 (atoms[i].mass>=14.0) && (atoms[i].mass<=18.0)){
3054 void Molecule::plgLoadBonds(
int *from,
int *to){
3055 bonds =
new Bond[numBonds];
3056 int realNumBonds = 0;
3057 for(
int i=0; i<numBonds; i++) {
3058 Bond *thisBond = bonds+realNumBonds;
3059 thisBond->
atom1 = from[i]-1;
3060 thisBond->
atom2 = to[i]-1;
3063 atomNames[thisBond->
atom1].atomtype,
3064 atomNames[thisBond->
atom2].atomtype,
3070 if (is_lonepairs_psf) {
3072 if(k!=0. || is_lp(thisBond->
atom1) ||
3073 is_lp(thisBond->
atom2)) {
3077 if(k != 0.) realNumBonds++;
3081 if(numBonds != realNumBonds) {
3082 iout <<
iWARN <<
"Ignored" << numBonds-realNumBonds <<
3083 "bonds with zero force constants.\n" <<
endi;
3084 iout <<
iWARN <<
"Will get H-H distance in rigid H20 from H-O-H angle.\n" <<
endi;
3086 numBonds = realNumBonds;
3089 void Molecule::plgLoadAngles(
int *plgAngles)
3091 angles=
new Angle[numAngles];
3092 int *atomid = plgAngles;
3093 int numRealAngles = 0;
3094 for(
int i=0; i<numAngles; i++) {
3095 Angle *thisAngle = angles+numRealAngles;
3096 thisAngle->
atom1 = atomid[0]-1;
3097 thisAngle->
atom2 = atomid[1]-1;
3098 thisAngle->
atom3 = atomid[2]-1;
3102 atomNames[thisAngle->
atom1].atomtype,
3103 atomNames[thisAngle->
atom2].atomtype,
3104 atomNames[thisAngle->
atom3].atomtype,
3107 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n" 3111 Real k, t0, k_ub, r_ub;
3112 if ( thisAngle->
angle_type == -1 ) { k = -1.; k_ub = -1.; }
else 3114 if(k!=0. || k_ub!=0.) numRealAngles++;
3117 if(numAngles != numRealAngles) {
3118 iout <<
iWARN <<
"Ignored" << numAngles-numRealAngles <<
3119 " angles with zero force constants.\n" <<
endi;
3121 numAngles = numRealAngles;
3124 void Molecule::plgLoadDihedrals(
int *plgDihedrals)
3126 std::map< std::string, int > cache;
3129 int multiplicity = 1;
3131 lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
3132 dihedrals =
new Dihedral[numDihedrals];
3133 int numRealDihedrals = 0;
3134 int *atomid = plgDihedrals;
3135 for(
int i=0; i<numDihedrals; i++, atomid+=4) {
3136 Dihedral *thisDihedral = dihedrals + numRealDihedrals;
3138 for(
int j=0; j<4; j++) {
3139 if(atomid[j] != lastAtomIds[j]) {
3140 duplicate_bond =
FALSE;
3142 lastAtomIds[j] = atomid[j];
3145 if(duplicate_bond) {
3147 if(multiplicity==2) {
3148 numMultipleDihedrals++;
3155 thisDihedral->
atom1 = atomid[0]-1;
3156 thisDihedral->
atom2 = atomid[1]-1;
3157 thisDihedral->
atom3 = atomid[2]-1;
3158 thisDihedral->
atom4 = atomid[3]-1;
3161 sprintf(query,
"%10s :: %10s :: %10s :: %10s :: %d",
3162 atomNames[atomid[0]-1].atomtype,
3163 atomNames[atomid[1]-1].atomtype,
3164 atomNames[atomid[2]-1].atomtype,
3165 atomNames[atomid[3]-1].atomtype,
3167 auto search = cache.find(query);
3168 if ( search != cache.end() ) {
3172 atomNames[atomid[0]-1].atomtype,
3173 atomNames[atomid[1]-1].atomtype,
3174 atomNames[atomid[2]-1].atomtype,
3175 atomNames[atomid[3]-1].atomtype,
3176 thisDihedral, multiplicity,
simParams->alchOn ? -1 : 0);
3178 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n" 3185 numDihedrals = numRealDihedrals;
3188 void Molecule::plgLoadImpropers(
int *plgImpropers)
3191 int multiplicity = 1;
3193 lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
3194 impropers =
new Improper[numImpropers];
3195 int numRealImpropers = 0;
3196 int *atomid = plgImpropers;
3197 for(
int i=0; i<numImpropers; i++, atomid+=4) {
3198 Improper *thisImproper = impropers + numRealImpropers;
3200 for(
int j=0; j<4; j++) {
3201 if(atomid[j] != lastAtomIds[j]) {
3202 duplicate_bond =
FALSE;
3204 lastAtomIds[j] = atomid[j];
3207 if(duplicate_bond) {
3209 if(multiplicity==2) {
3210 numMultipleImpropers++;
3217 thisImproper->
atom1 = atomid[0]-1;
3218 thisImproper->
atom2 = atomid[1]-1;
3219 thisImproper->
atom3 = atomid[2]-1;
3220 thisImproper->
atom4 = atomid[3]-1;
3223 atomNames[atomid[0]-1].atomtype,
3224 atomNames[atomid[1]-1].atomtype,
3225 atomNames[atomid[2]-1].atomtype,
3226 atomNames[atomid[3]-1].atomtype,
3227 thisImproper, multiplicity);
3230 numImpropers = numRealImpropers;
3233 void Molecule::plgLoadCrossterms(
int *plgCterms)
3237 for(
int i=0; i<8; i++)
3240 crossterms =
new Crossterm[numCrossterms];
3241 int numRealCrossterms = 0;
3242 int *atomid = plgCterms;
3243 for(
int i=0; i<numCrossterms; i++, atomid+=8) {
3244 Crossterm *thisCrossterm = crossterms + numRealCrossterms;
3246 for(
int j=0; j<8; j++) {
3247 if(atomid[j] != lastAtomIds[j]) {
3248 duplicate_bond =
FALSE;
3250 lastAtomIds[j] = atomid[j];
3253 if(duplicate_bond) {
3256 numRealCrossterms++;
3258 thisCrossterm->
atom1 = atomid[0]-1;
3259 thisCrossterm->
atom2 = atomid[1]-1;
3260 thisCrossterm->
atom3 = atomid[2]-1;
3261 thisCrossterm->
atom4 = atomid[3]-1;
3262 thisCrossterm->
atom5 = atomid[4]-1;
3263 thisCrossterm->
atom6 = atomid[5]-1;
3264 thisCrossterm->
atom7 = atomid[6]-1;
3265 thisCrossterm->
atom8 = atomid[7]-1;
3268 atomNames[atomid[0]-1].atomtype,
3269 atomNames[atomid[1]-1].atomtype,
3270 atomNames[atomid[2]-1].atomtype,
3271 atomNames[atomid[3]-1].atomtype,
3272 atomNames[atomid[4]-1].atomtype,
3273 atomNames[atomid[5]-1].atomtype,
3274 atomNames[atomid[6]-1].atomtype,
3275 atomNames[atomid[7]-1].atomtype,
3279 numCrossterms = numRealCrossterms;
3283 occupancy =
new float[numAtoms];
3284 for(
int i=0; i<numAtoms; i++) {
3285 occupancy[i] = atomarray[i].occupancy;
3290 bfactor =
new float[numAtoms];
3291 for(
int i=0; i<numAtoms; i++) {
3292 bfactor[i] = atomarray[i].bfactor;
3300 return v1.size() > v2.size();
3310 vector<vector<int> > atomBonds(numAtoms);
3311 for (i = 0; i < numRealBonds; i++) {
3312 int a = bonds[i].atom1;
3313 int b = bonds[i].atom2;
3314 atomBonds[a].push_back(b);
3315 atomBonds[b].push_back(a);
3318 vector<int> atomsInMolecules(numAtoms, -1);
3320 for (i = 0; i < numAtoms; i++) {
3321 if (atomsInMolecules[i] == -1) {
3322 vector<int> atomList, neighborList;
3323 atomList.push_back(i);
3324 neighborList.push_back(0);
3325 int currentMolecule = molNumber++;
3327 while (atomList.size() > 0) {
3328 int particle = atomList.back();
3329 atomsInMolecules[particle] = currentMolecule;
3330 int& neighbor = neighborList.back();
3331 while (neighbor < atomBonds[particle].size() && atomsInMolecules[atomBonds[particle][neighbor]] != -1) {
3335 if (neighbor < atomBonds[particle].size()) {
3336 atomList.push_back(atomBonds[particle][neighbor]);
3337 neighborList.push_back(0);
3339 atomList.pop_back();
3340 neighborList.pop_back();
3346 numMolecules = molNumber;
3348 vector<vector<int> > molecules(numMolecules);
3349 for (i = 0; i < numAtoms; i++) {
3350 molecules[atomsInMolecules[i]].push_back(i);
3355 sort(molecules.begin(), molecules.end(),
sizeColumn);
3358 numLargeMolecules = 0;
3359 for (i = 0; i < numMolecules; i++){
3361 ++numLargeMolecules;
3373 moleculeStartIndex =
new int32[numMolecules + 1];
3374 moleculeAtom =
new int32[numAtoms];
3376 for (i = 0; i < numMolecules; i++) {
3377 moleculeStartIndex[i] = index;
3378 for (j = 0; j < molecules[i].size(); j++) {
3379 moleculeAtom[index++] = molecules[i][j];
3383 moleculeStartIndex[numMolecules] = index;
3386 printf(
"Number of Molecule with atom greater than %5d: %5d \n: ",
LARGEMOLTH, numLargeMolecules);
3387 for (i = 0; i < numMolecules; i++) {
3388 int startIdx = moleculeStartIndex[i];
3389 int endIdx = moleculeStartIndex[i+1];
3390 printf(
"Molecule Idx %5d: ", i);
3391 for (j = startIdx; j < endIdx; j++) {
3392 int atom = moleculeAtom[j];
3393 printf(
"%7d ", atom);
3411 void Molecule::build_lists_by_atom()
3415 register int numFixedAtoms = this->numFixedAtoms;
3418 if (
simParams->fixedAtomsForces ) numFixedAtoms = 0;
3425 bondsWithAtom =
new int32 *[numAtoms];
3426 cluster =
new int32 [numAtoms];
3427 clusterSize =
new int32 [numAtoms];
3429 bondsByAtom =
new int32 *[numAtoms];
3430 anglesByAtom =
new int32 *[numAtoms];
3431 dihedralsByAtom =
new int32 *[numAtoms];
3432 impropersByAtom =
new int32 *[numAtoms];
3433 crosstermsByAtom =
new int32 *[numAtoms];
3435 exclusionsByAtom =
new int32 *[numAtoms];
3436 fullExclusionsByAtom =
new int32 *[numAtoms];
3437 modExclusionsByAtom =
new int32 *[numAtoms];
3440 gromacsPairByAtom =
new int32 *[numAtoms];
3445 const int pair_self =
3448 DebugM(3,
"Building bond lists.\n");
3451 for (i=0; i<numAtoms; i++)
3455 for (i=0; i<numRealBonds; i++)
3457 byAtomSize[bonds[i].atom1]++;
3458 byAtomSize[bonds[i].atom2]++;
3460 for (i=0; i<numAtoms; i++)
3462 bondsWithAtom[i] = tmpArena->getNewArray(byAtomSize[i]+1);
3463 bondsWithAtom[i][byAtomSize[i]] = -1;
3466 for (i=0; i<numRealBonds; i++)
3468 int a1 = bonds[i].atom1;
3469 int a2 = bonds[i].atom2;
3470 bondsWithAtom[a1][byAtomSize[a1]++] = i;
3471 bondsWithAtom[a2][byAtomSize[a2]++] = i;
3479 DebugM(3,
"Calculating exclusions for QM simulation.\n");
3482 delete_qm_bonded() ;
3484 DebugM(3,
"Re-Building bond lists.\n");
3488 for (i=0; i<numAtoms; i++)
3492 for (i=0; i<numRealBonds; i++)
3494 byAtomSize[bonds[i].atom1]++;
3495 byAtomSize[bonds[i].atom2]++;
3497 for (i=0; i<numAtoms; i++)
3499 bondsWithAtom[i][byAtomSize[i]] = -1;
3502 for (i=0; i<numRealBonds; i++)
3504 int a1 = bonds[i].atom1;
3505 int a2 = bonds[i].atom2;
3506 bondsWithAtom[a1][byAtomSize[a1]++] = i;
3507 bondsWithAtom[a2][byAtomSize[a2]++] = i;
3512 for (i=0; i<numAtoms; i++) {
3515 for (i=0; i<numAtoms; i++) {
3517 while ( cluster[ci] != ci ) ci = cluster[ci];
3518 for (
int32 *b = bondsWithAtom[i]; *b != -1; ++b ) {
3519 int a = bonds[*b].atom1;
3520 if ( a == i ) a = bonds[*b].atom2;
3523 while ( cluster[ca] != ca ) ca = cluster[ca];
3524 if ( ca > ci ) cluster[ca] = cluster[ci];
3525 else cluster[ci] = cluster[ca];
3531 for (i=0; i<numAtoms; i++) {
3532 int ci = cluster[i];
3533 if ( cluster[ci] != ci ) {
3535 cluster[i] = cluster[ci];
3541 for (i=0; i<numAtoms; i++) {
3544 for (i=0; i<numAtoms; i++) {
3545 clusterSize[cluster[i]] += 1;
3558 for (i=0; i<numAtoms; i++)
3563 for (i=0; i<numBonds; i++)
3565 if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3566 && fixedAtomFlags[bonds[i].atom2] )
continue;
3568 if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1)
continue;
3569 byAtomSize[bonds[i].atom1]++;
3572 for (i=0; i<numAtoms; i++)
3574 bondsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3575 bondsByAtom[i][byAtomSize[i]] = -1;
3578 for (i=0; i<numBonds; i++)
3580 if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3581 && fixedAtomFlags[bonds[i].atom2] )
continue;
3582 if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1)
continue;
3583 int a1 = bonds[i].atom1;
3584 bondsByAtom[a1][byAtomSize[a1]++] = i;
3586 for (i=0; i<numBonds; ++i) {
3587 int a1 = bonds[i].atom1;
3588 int a2 = bonds[i].atom2;
3592 sprintf(buff,
"Atom %d is bonded to itself", a1+1);
3595 for ( j = 0; j < byAtomSize[a1]; ++j ) {
3596 int b = bondsByAtom[a1][j];
3597 int ba1 = bonds[b].atom1;
3598 int ba2 = bonds[b].atom2;
3599 if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3601 sprintf(buff,
"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3605 for ( j = 0; j < byAtomSize[a2]; ++j ) {
3606 int b = bondsByAtom[a2][j];
3607 int ba1 = bonds[b].atom1;
3608 int ba2 = bonds[b].atom2;
3609 if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3611 sprintf(buff,
"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3617 DebugM(3,
"Building angle lists.\n");
3620 for (i=0; i<numAtoms; i++)
3625 for (i=0; i<numAngles; i++)
3627 if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3628 && fixedAtomFlags[angles[i].atom2]
3629 && fixedAtomFlags[angles[i].atom3] )
continue;
3630 if ( pair_self && fepAtomFlags[angles[i].atom1] != 1)
continue;
3631 byAtomSize[angles[i].atom1]++;
3634 for (i=0; i<numAtoms; i++)
3636 anglesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3637 anglesByAtom[i][byAtomSize[i]] = -1;
3640 for (i=0; i<numAngles; i++)
3642 if ( pair_self && fepAtomFlags[angles[i].atom1] != 1)
continue;
3643 if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3644 && fixedAtomFlags[angles[i].atom2]
3645 && fixedAtomFlags[angles[i].atom3] )
continue;
3646 int a1 = angles[i].atom1;
3647 anglesByAtom[a1][byAtomSize[a1]++] = i;
3650 DebugM(3,
"Building improper lists.\n");
3653 for (i=0; i<numAtoms; i++)
3657 numCalcImpropers = 0;
3658 for (i=0; i<numImpropers; i++)
3660 if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3661 && fixedAtomFlags[impropers[i].atom2]
3662 && fixedAtomFlags[impropers[i].atom3]
3663 && fixedAtomFlags[impropers[i].atom4] )
continue;
3664 if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1)
continue;
3665 byAtomSize[impropers[i].atom1]++;
3668 for (i=0; i<numAtoms; i++)
3670 impropersByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3671 impropersByAtom[i][byAtomSize[i]] = -1;
3674 for (i=0; i<numImpropers; i++)
3676 if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3677 && fixedAtomFlags[impropers[i].atom2]
3678 && fixedAtomFlags[impropers[i].atom3]
3679 && fixedAtomFlags[impropers[i].atom4] )
continue;
3680 if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1)
continue;
3681 int a1 = impropers[i].atom1;
3682 impropersByAtom[a1][byAtomSize[a1]++] = i;
3685 DebugM(3,
"Building dihedral lists.\n");
3688 for (i=0; i<numAtoms; i++)
3692 numCalcDihedrals = 0;
3693 for (i=0; i<numDihedrals; i++)
3695 if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3696 && fixedAtomFlags[dihedrals[i].atom2]
3697 && fixedAtomFlags[dihedrals[i].atom3]
3698 && fixedAtomFlags[dihedrals[i].atom4] )
continue;
3699 if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1)
continue;
3700 byAtomSize[dihedrals[i].atom1]++;
3703 for (i=0; i<numAtoms; i++)
3705 dihedralsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3706 dihedralsByAtom[i][byAtomSize[i]] = -1;
3709 for (i=0; i<numDihedrals; i++)
3711 if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3712 && fixedAtomFlags[dihedrals[i].atom2]
3713 && fixedAtomFlags[dihedrals[i].atom3]
3714 && fixedAtomFlags[dihedrals[i].atom4] )
continue;
3715 if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1)
continue;
3716 int a1 = dihedrals[i].atom1;
3717 dihedralsByAtom[a1][byAtomSize[a1]++] = i;
3720 DebugM(3,
"Building crossterm lists.\n");
3723 for (i=0; i<numAtoms; i++)
3727 numCalcCrossterms = 0;
3728 for (i=0; i<numCrossterms; i++)
3730 if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3731 && fixedAtomFlags[crossterms[i].atom2]
3732 && fixedAtomFlags[crossterms[i].atom3]
3733 && fixedAtomFlags[crossterms[i].atom4]
3734 && fixedAtomFlags[crossterms[i].atom5]
3735 && fixedAtomFlags[crossterms[i].atom6]
3736 && fixedAtomFlags[crossterms[i].atom7]
3737 && fixedAtomFlags[crossterms[i].atom8] )
continue;
3738 if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1)
continue;
3739 byAtomSize[crossterms[i].atom1]++;
3740 numCalcCrossterms++;
3742 for (i=0; i<numAtoms; i++)
3744 crosstermsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3745 crosstermsByAtom[i][byAtomSize[i]] = -1;
3748 for (i=0; i<numCrossterms; i++)
3750 if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3751 && fixedAtomFlags[crossterms[i].atom2]
3752 && fixedAtomFlags[crossterms[i].atom3]
3753 && fixedAtomFlags[crossterms[i].atom4]
3754 && fixedAtomFlags[crossterms[i].atom5]
3755 && fixedAtomFlags[crossterms[i].atom6]
3756 && fixedAtomFlags[crossterms[i].atom7]
3757 && fixedAtomFlags[crossterms[i].atom8] )
continue;
3758 if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1)
continue;
3759 int a1 = crossterms[i].atom1;
3760 crosstermsByAtom[a1][byAtomSize[a1]++] = i;
3764 if (is_lonepairs_psf) {
3766 DebugM(3,
"Initializing lone pair host index array.\n");
3767 lphostIndexes =
new int32[numAtoms];
3768 for (i = 0; i < numAtoms; i++) {
3769 lphostIndexes[i] = -1;
3771 for (i = 0; i < numLphosts; i++) {
3772 int32 index = lphosts[i].atom1;
3773 lphostIndexes[index] = i;
3779 DebugM(3,
"Building gromacsPair lists.\n");
3782 for (i=0; i<numAtoms; i++)
3789 if ( numFixedAtoms && fixedAtomFlags[
gromacsPair[i].atom1]
3790 && fixedAtomFlags[
gromacsPair[i].atom2] )
continue;
3791 if ( pair_self && fepAtomFlags[
gromacsPair[i].atom1] != 1)
continue;
3795 for (i=0; i<numAtoms; i++)
3797 gromacsPairByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3798 gromacsPairByAtom[i][byAtomSize[i]] = -1;
3803 if ( numFixedAtoms && fixedAtomFlags[
gromacsPair[i].atom1]
3804 && fixedAtomFlags[
gromacsPair[i].atom2] )
continue;
3805 if ( pair_self && fepAtomFlags[
gromacsPair[i].atom1] != 1)
continue;
3807 gromacsPairByAtom[a1][byAtomSize[a1]++] = i;
3812 DebugM(3,
"Building exclusion data.\n");
3819 delete [] bondsWithAtom; bondsWithAtom = 0;
3820 delete tmpArena; tmpArena = 0;
3822 if (exclusions != NULL)
3823 delete [] exclusions;
3826 numTotalExclusions = exclusionSet.size();
3828 iout <<
iINFO <<
"ADDED " << (numTotalExclusions - numExclusions) <<
" IMPLICIT EXCLUSIONS\n" <<
endi;
3830 exclusions =
new Exclusion[numTotalExclusions];
3832 for ( exclIter=exclIter.begin(),i=0; exclIter != exclIter.end(); exclIter++,i++ )
3834 exclusions[i] = *exclIter;
3838 exclusionSet.clear();
3840 DebugM(3,
"Building exclusion lists.\n");
3842 for (i=0; i<numAtoms; i++)
3846 numCalcExclusions = 0;
3847 numCalcFullExclusions = 0;
3848 for (i=0; i<numTotalExclusions; i++)
3850 if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3851 && fixedAtomFlags[exclusions[i].atom2] )
continue;
3852 byAtomSize[exclusions[i].atom1]++;
3853 numCalcExclusions++;
3854 if ( ! exclusions[i].modified ) numCalcFullExclusions++;
3857 for (i=0; i<numAtoms; i++)
3859 exclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3860 exclusionsByAtom[i][byAtomSize[i]] = -1;
3863 for (i=0; i<numTotalExclusions; i++)
3865 if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3866 && fixedAtomFlags[exclusions[i].atom2] )
continue;
3867 int a1 = exclusions[i].atom1;
3868 exclusionsByAtom[a1][byAtomSize[a1]++] = i;
3873 for (i=0; i<numAtoms; i++)
3879 for (i=0; i<numTotalExclusions; i++)
3881 if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3882 && fixedAtomFlags[exclusions[i].atom2] )
continue;
3883 if ( exclusions[i].modified ) {
3884 byAtomSize2[exclusions[i].atom1]++;
3885 byAtomSize2[exclusions[i].atom2]++;
3887 byAtomSize[exclusions[i].atom1]++;
3888 byAtomSize[exclusions[i].atom2]++;
3892 for (i=0; i<numAtoms; i++)
3894 fullExclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3895 fullExclusionsByAtom[i][0] = 0;
3896 modExclusionsByAtom[i] = arena->getNewArray(byAtomSize2[i]+1);
3897 modExclusionsByAtom[i][0] = 0;
3900 for (i=0; i<numTotalExclusions; i++)
3902 int a1 = exclusions[i].atom1;
3903 int a2 = exclusions[i].atom2;
3904 if ( numFixedAtoms && fixedAtomFlags[a1]
3905 && fixedAtomFlags[a2] )
continue;
3907 if ( exclusions[i].modified ) {
3908 l1 = modExclusionsByAtom[a1];
3909 l2 = modExclusionsByAtom[a2];
3911 l1 = fullExclusionsByAtom[a1];
3912 l2 = fullExclusionsByAtom[a2];
3918 if ( ! CkMyPe() &&
simParams->printExclusions ) {
3919 for (i=0; i<numAtoms; i++) {
3920 int32 *lf = fullExclusionsByAtom[i];
3921 iout <<
"EXCL " << i <<
" FULL";
3923 for (
int j = 0; j < nf; ++j ) {
3924 iout <<
" " << *(lf++);
3927 int32 *lm = modExclusionsByAtom[i];
3928 iout <<
"EXCL " << i <<
" MOD";
3930 for (
int j = 0; j < nm; ++j ) {
3931 iout <<
" " << *(lm++);
3938 if (is_drude_psf ||
simParams->drudeOn) {
3944 if (tholes != NULL)
delete[] tholes;
3947 if (oneFourNbTholes != NULL)
delete[] oneFourNbTholes;
3948 numOneFourNbTholes = 0;
3951 for (i = 0; i < numTotalExclusions; i++) {
3953 if (exclusions[i].modified)
continue;
3954 int a1 = exclusions[i].atom1;
3955 int a2 = exclusions[i].atom2;
3956 if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3962 if (numTholes != 0) tholes =
new Thole[numTholes];
3969 std::vector<OneFourNbThole> oneFourNbTholesTmp;
3971 for (i = 0; i < numTotalExclusions; i++) {
3972 int a1 = exclusions[i].atom1;
3973 int a2 = exclusions[i].atom2;
3974 if (exclusions[i].modified) {
3975 if (a1 < numAtoms-1 && a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3979 if (((nbthole_array[k].ind1 == atoms[a1].vdw_type) && (nbthole_array[k].ind2 == atoms[a2].vdw_type)) ||
3980 ((nbthole_array[k].
ind2 == atoms[a1].vdw_type) && (nbthole_array[k].ind1 == atoms[a2].vdw_type))) {
3982 const Real aprod = drudeConsts[a1].alpha * drudeConsts[a2].alpha;
3983 const Real apower = (aprod <= 0 ? 0 : powf(aprod, -1.f/6));
3984 aa = apower * nbthole_array[k].
tholeij;
3989 oneFourNbTholesTmp.push_back(
OneFourNbThole{a1, a1+1, a2, a2+1, aa});
3990 DebugM(3,
"Add 1-4 NbThole term for (" << a1 <<
", " << a1+1 <<
", " << a2 <<
", " << a2+1 <<
")\n");
3995 if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3996 Real thsum = drudeConsts[a1].thole + drudeConsts[a2].thole;
3997 Real aprod = drudeConsts[a1].alpha * drudeConsts[a2].alpha;
3999 Real apower = (aprod <= 0 ? 0 : powf(aprod, -1.f/6));
4000 tholes[nt].atom1 = a1;
4001 tholes[nt].atom2 = a1+1;
4002 tholes[nt].atom3 = a2;
4003 tholes[nt].atom4 = a2+1;
4004 tholes[nt].aa = apower * thsum;
4005 tholes[nt].qq = c * atoms[a1+1].charge * atoms[a2+1].charge;
4010 numOneFourNbTholes = oneFourNbTholesTmp.size();
4011 if (!oneFourNbTholesTmp.empty()) {
4013 std::copy(oneFourNbTholesTmp.begin(), oneFourNbTholesTmp.end(), oneFourNbTholes);
4017 DebugM(3,
"Building Thole correction term lists.\n");
4018 tholesByAtom =
new int32 *[numAtoms];
4020 for (i = 0; i < numAtoms; i++) {
4024 for (i = 0; i < numTholes; i++) {
4025 if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
4026 && fixedAtomFlags[tholes[i].atom2]
4027 && fixedAtomFlags[tholes[i].atom3]
4028 && fixedAtomFlags[tholes[i].atom4] )
continue;
4029 if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1)
continue;
4030 byAtomSize[tholes[i].atom1]++;
4033 for (i = 0; i < numAtoms; i++) {
4034 tholesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
4035 tholesByAtom[i][byAtomSize[i]] = -1;
4038 for (i = 0; i < numTholes; i++) {
4039 if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
4040 && fixedAtomFlags[tholes[i].atom2]
4041 && fixedAtomFlags[tholes[i].atom3]
4042 && fixedAtomFlags[tholes[i].atom4] )
continue;
4043 if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1)
continue;
4044 int a1 = tholes[i].atom1;
4045 tholesByAtom[a1][byAtomSize[a1]++] = i;
4049 oneFourNbTholesByAtom =
new int32 *[numAtoms];
4050 for (i = 0; i <numAtoms; ++i) {
4053 numCalcOneFourNbTholes = 0;
4054 for (i = 0; i < numOneFourNbTholes; ++i) {
4055 if ( numFixedAtoms && fixedAtomFlags[oneFourNbTholes[i].atom1]
4056 && fixedAtomFlags[oneFourNbTholes[i].atom2]
4057 && fixedAtomFlags[oneFourNbTholes[i].atom3]
4058 && fixedAtomFlags[oneFourNbTholes[i].atom4] )
continue;
4059 if ( pair_self && fepAtomFlags[oneFourNbTholes[i].atom1] != 1)
continue;
4060 byAtomSize[oneFourNbTholes[i].atom1]++;
4061 numCalcOneFourNbTholes++;
4063 for (i = 0; i < numAtoms; i++) {
4064 oneFourNbTholesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
4065 oneFourNbTholesByAtom[i][byAtomSize[i]] = -1;
4068 for (i = 0; i < numOneFourNbTholes; i++) {
4069 if ( numFixedAtoms && fixedAtomFlags[oneFourNbTholes[i].atom1]
4070 && fixedAtomFlags[oneFourNbTholes[i].atom2]
4071 && fixedAtomFlags[oneFourNbTholes[i].atom3]
4072 && fixedAtomFlags[oneFourNbTholes[i].atom4] )
continue;
4073 if ( pair_self && fepAtomFlags[oneFourNbTholes[i].atom1] != 1)
continue;
4074 int a1 = oneFourNbTholes[i].atom1;
4075 oneFourNbTholesByAtom[a1][byAtomSize[a1]++] = i;
4079 DebugM(3,
"Building anisotropic term lists.\n");
4080 anisosByAtom =
new int32 *[numAtoms];
4082 for (i = 0; i < numAtoms; i++) {
4086 for (i = 0; i < numAnisos; i++) {
4087 if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
4088 && fixedAtomFlags[anisos[i].atom2]
4089 && fixedAtomFlags[anisos[i].atom3]
4090 && fixedAtomFlags[anisos[i].atom4] )
continue;
4091 if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1)
continue;
4092 byAtomSize[anisos[i].atom1]++;
4095 for (i = 0; i < numAtoms; i++) {
4096 anisosByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
4097 anisosByAtom[i][byAtomSize[i]] = -1;
4100 for (i = 0; i < numAnisos; i++) {
4101 if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
4102 && fixedAtomFlags[anisos[i].atom2]
4103 && fixedAtomFlags[anisos[i].atom3]
4104 && fixedAtomFlags[anisos[i].atom4] )
continue;
4105 if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1)
continue;
4106 int a1 = anisos[i].atom1;
4107 anisosByAtom[a1][byAtomSize[a1]++] = i;
4109 iout <<
iINFO <<
"There are " << numTholes <<
" Thole terms.\n" <<
endi;
4110 iout <<
iINFO <<
"There are " << numAnisos <<
" Anisotropic polarization terms.\n" <<
endi;
4111 iout <<
iINFO <<
"There are " << numOneFourNbTholes <<
" 1-4 NbThole terms.\n" <<
endi;
4115 delete [] byAtomSize; byAtomSize = 0;
4116 delete [] byAtomSize2; byAtomSize2 = 0;
4122 for (i=0; i<numAtoms; i++)
4124 all_exclusions[i].
min = numAtoms;
4125 all_exclusions[i].max = -1;
4127 for (i=0; i<numTotalExclusions; i++)
4130 int a1 = exclusions[i].atom1;
4131 int a2 = exclusions[i].atom2;
4132 if ( numFixedAtoms && fixedAtomFlags[a1]
4133 && fixedAtomFlags[a2] )
continue;
4134 if ( all_exclusions[a1].min > a2 ) all_exclusions[a1].min = a2;
4135 if ( all_exclusions[a2].min > a1 ) all_exclusions[a2].min = a1;
4136 if ( a2 > all_exclusions[a1].max ) all_exclusions[a1].max = a2;
4137 if ( a1 > all_exclusions[a2].max ) all_exclusions[a2].max = a1;
4141 int maxDenseAllFull = 0;
4142 int numDenseAllFull = 0;
4143 for (i=0; i<numAtoms; i++) {
4144 int iInMiddle = ( i < all_exclusions[i].max &&
4145 i > all_exclusions[i].min ) ? 1 : 0;
4146 int s = all_exclusions[i].max - all_exclusions[i].min + 1;
4147 if ( s == fullExclusionsByAtom[i][0] + iInMiddle ) {
4148 if ( s > maxDenseAllFull ) maxDenseAllFull = s;
4149 all_exclusions[i].flags = (
char*)-1;
4151 all_exclusions[i].flags = 0;
4154 char *denseFullArray = exclArena->getNewArray(maxDenseAllFull);
4155 for ( i=0; i<maxDenseAllFull; ++i ) denseFullArray[i] =
EXCHCK_FULL;
4157 int exclmem = maxDenseAllFull;
4158 int maxExclusionFlags =
simParams->maxExclusionFlags;
4159 for (i=0; i<numAtoms; i++) {
4160 int s = all_exclusions[i].max - all_exclusions[i].min + 1;
4161 if ( all_exclusions[i].max != -1 ) {
4162 if ( all_exclusions[i].flags ) {
4163 all_exclusions[i].flags = denseFullArray;
4165 }
else if ( s < maxExclusionFlags ) {
4166 char *f = all_exclusions[i].flags = exclArena->getNewArray(s);
4167 for (
int k=0; k<s; ++k ) f[k] = 0;
4170 all_exclusions[i].flags = 0;
4173 all_exclusions[i].flags = (
char*)-1;
4177 iout <<
iINFO << numTotalExclusions <<
" exclusions consume " 4178 << exclmem <<
" bytes.\n" <<
endi;
4180 <<
" atoms sharing one array.\n" <<
endi;
4182 for (i=0; i<numTotalExclusions; i++)
4184 int a1 = exclusions[i].atom1;
4185 int a2 = exclusions[i].atom2;
4186 if ( numFixedAtoms && fixedAtomFlags[a1]
4187 && fixedAtomFlags[a2] )
continue;
4188 if ( exclusions[i].modified ) {
4189 if ( all_exclusions[a1].flags )
4190 all_exclusions[a1].flags[a2-all_exclusions[a1].min] =
EXCHCK_MOD;
4191 if ( all_exclusions[a2].flags )
4192 all_exclusions[a2].flags[a1-all_exclusions[a2].min] =
EXCHCK_MOD;
4194 if ( all_exclusions[a1].flags )
4195 all_exclusions[a1].flags[a2-all_exclusions[a1].min] =
EXCHCK_FULL;
4196 if ( all_exclusions[a2].flags )
4197 all_exclusions[a2].flags[a1-all_exclusions[a2].min] =
EXCHCK_FULL;
4226 void Molecule::build_exclusions()
4234 for (i=0; i<numExclusions; i++)
4236 exclusionSet.add(exclusions[i]);
4244 switch (exclude_flag)
4271 if (is_lonepairs_psf || is_drude_psf) {
4272 build_inherited_excl(
SCALED14 == exclude_flag);
4291 void Molecule::build_inherited_excl(
int modified) {
4293 int32 *bond1, *bond2, *bond3, *bond4, *bond5;
4294 int32 i, j, mid1, mid2, mid3, mid4;
4298 for (i = 0; i < numAtoms; i++) {
4300 if (!is_drude(i) && !is_lp(i))
continue;
4304 bond1 = bondsWithAtom[i];
4308 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4309 sprintf(err_msg,
"FOUND ISOLATED %s PARTICLE %d", idescrip, i+1);
4312 if (-1 != *(bond1+1)) {
4314 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4315 sprintf(err_msg,
"FOUND MULTIPLY LINKED %s PARTICLE %d",
4321 mid1 = bonds[*bond1].atom1;
4322 if (mid1 == i) mid1 = bonds[*bond1].atom2;
4325 if (is_drude(mid1) || is_lp(mid1) || is_hydrogen(mid1)) {
4327 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4328 sprintf(err_msg,
"PARENT ATOM %d of %s PARTICLE %d " 4329 "IS NOT HEAVY ATOM", mid1+1, idescrip, i+1);
4333 if (exclude_flag ==
NONE) {
4343 bond2 = bondsWithAtom[mid1];
4344 while (*bond2 != -1) {
4345 j = bonds[*bond2].atom1;
4346 if ((is_drude(j) || is_lp(j)) && j != mid1) {
4347 if (i < j) exclusionSet.add(
Exclusion(i, j));
4348 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4350 j = bonds[*bond2].atom2;
4351 if ((is_drude(j) || is_lp(j)) && j != mid1) {
4352 if (i < j) exclusionSet.add(
Exclusion(i, j));
4353 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4361 bond2 = bondsWithAtom[mid1];
4364 while (*bond2 != -1) {
4365 if (bonds[*bond2].atom1 == mid1) {
4366 mid2 = bonds[*bond2].atom2;
4369 mid2 = bonds[*bond2].atom1;
4379 if (exclude_flag ==
ONETWO) {
4389 bond3 = bondsWithAtom[mid2];
4390 while (*bond3 != -1) {
4391 j = bonds[*bond3].atom1;
4392 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4393 if (i < j) exclusionSet.add(
Exclusion(i, j));
4394 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4396 j = bonds[*bond3].atom2;
4397 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4398 if (i < j) exclusionSet.add(
Exclusion(i, j));
4399 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4407 bond3 = bondsWithAtom[mid2];
4410 while (*bond3 != -1) {
4412 if (bonds[*bond3].atom1 == mid2) {
4413 mid3 = bonds[*bond3].atom2;
4416 mid3 = bonds[*bond3].atom1;
4430 else if (mid3 < i) {
4436 bond4 = bondsWithAtom[mid3];
4437 while (*bond4 != -1) {
4438 j = bonds[*bond4].atom1;
4439 if ((is_drude(j) || is_lp(j)) && j != mid3) {
4440 if (i < j) exclusionSet.add(
Exclusion(i, j));
4441 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4443 j = bonds[*bond4].atom2;
4444 if ((is_drude(j) || is_lp(j)) && j != mid3) {
4445 if (i < j) exclusionSet.add(
Exclusion(i, j));
4446 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4454 bond4 = bondsWithAtom[mid3];
4457 while (*bond4 != -1) {
4459 if (bonds[*bond4].atom1 == mid3) {
4460 mid4 = bonds[*bond4].atom2;
4463 mid4 = bonds[*bond4].atom1;
4473 if (is_drude(mid4) || is_lp(mid4)) {
4478 else if (mid4 < i) {
4489 int modi = modified;
4491 int amin = (mid1 < mid4 ? mid1 : mid4);
4492 int amax = (mid1 >= mid4 ? mid1 : mid4);
4504 exclusionSet.add(
Exclusion(i, mid4, modi));
4506 else if (mid4 < i) {
4507 exclusionSet.add(
Exclusion(mid4, i, modi));
4512 bond5 = bondsWithAtom[mid4];
4513 while (*bond5 != -1) {
4514 j = bonds[*bond5].atom1;
4515 if ((is_drude(j) || is_lp(j)) && j != mid4) {
4516 if (i<j) exclusionSet.add(
Exclusion(i,j,modi));
4517 else if (j<i) exclusionSet.add(
Exclusion(j,i,modi));
4519 j = bonds[*bond5].atom2;
4520 if ((is_drude(j) || is_lp(j)) && j != mid4) {
4521 if (i<j) exclusionSet.add(
Exclusion(i,j,modi));
4522 else if (j<i) exclusionSet.add(
Exclusion(j,i,modi));
4552 void Molecule::build12excl(
void)
4558 for (i=0; i<numAtoms; i++)
4560 current_val = bondsWithAtom[i];
4563 while (*current_val != -1)
4565 if (bonds[*current_val].atom1 == i)
4567 if (i<bonds[*current_val].atom2)
4569 exclusionSet.add(
Exclusion(i,bonds[*current_val].atom2));
4574 if (i<bonds[*current_val].atom1)
4576 exclusionSet.add(
Exclusion(i,bonds[*current_val].atom1));
4592 void Molecule::build13excl(
void)
4594 int32 *bond1, *bond2;
4600 for (i=0; i<numAtoms; i++)
4602 bond1 = bondsWithAtom[i];
4605 while (*bond1 != -1)
4607 if (bonds[*bond1].atom1 == i)
4609 middle_atom=bonds[*bond1].atom2;
4613 middle_atom=bonds[*bond1].atom1;
4616 bond2 = bondsWithAtom[middle_atom];
4620 while (*bond2 != -1)
4622 if (bonds[*bond2].atom1 == middle_atom)
4624 if (i < bonds[*bond2].atom2)
4626 exclusionSet.add(
Exclusion(i,bonds[*bond2].atom2));
4631 if (i < bonds[*bond2].atom1)
4633 exclusionSet.add(
Exclusion(i,bonds[*bond2].atom1));
4653 void Molecule::build14excl(
int modified)
4655 int32 *bond1, *bond2, *bond3;
4660 for (i=0; i<numAtoms; i++)
4662 if (is_drude(i) || is_lp(i))
continue;
4665 bond1 = bondsWithAtom[i];
4667 while (*bond1 != -1)
4669 if (bonds[*bond1].atom1 == i)
4671 mid1=bonds[*bond1].atom2;
4675 mid1=bonds[*bond1].atom1;
4678 bond2 = bondsWithAtom[mid1];
4681 while (*bond2 != -1)
4683 if (bonds[*bond2].atom1 == mid1)
4685 mid2 = bonds[*bond2].atom2;
4689 mid2 = bonds[*bond2].atom1;
4701 bond3=bondsWithAtom[mid2];
4704 while (*bond3 != -1)
4706 if (bonds[*bond3].atom1 == mid2)
4708 int j = bonds[*bond3].atom2;
4714 if (i < j && !is_drude(j) && !is_lp(j))
4716 exclusionSet.add(
Exclusion(i,j,modified));
4721 int j = bonds[*bond3].atom1;
4727 if (i < j && !is_drude(j) && !is_lp(j))
4729 exclusionSet.add(
Exclusion(i,j,modified));
4751 void Molecule::stripFepExcl(
void)
4757 for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
4759 int t1 = get_fep_type(exclIter->atom1);
4760 int t2 = get_fep_type(exclIter->atom2);
4761 if (t1 && t2 && t1 !=t2 && abs(t1-t2) != 2) {
4762 fepExclusionSet.
add(*exclIter);
4765 }
else if (
simParams->pairInteractionOn ) {
4766 for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
4768 int ifep_type = get_fep_type(exclIter->atom1);
4769 int jfep_type = get_fep_type(exclIter->atom2);
4772 if (ifep_type != 1 || jfep_type != 1) {
4773 fepExclusionSet.
add(*exclIter);
4778 if (!(ifep_type == 1 && jfep_type == 2) &&
4779 !(ifep_type == 2 && jfep_type == 1)) {
4780 fepExclusionSet.
add(*exclIter);
4787 for ( fepIter=fepIter.begin(); fepIter != fepIter.end(); fepIter++ )
4789 exclusionSet.del(*fepIter);
4803 if ( (psf_file =
Fopen(fname,
"r")) == NULL)
4806 sprintf(err_msg,
"UNABLE TO OPEN THE COMPRESSED .psf FILE %s", fname);
4815 NAMD_die(
"The compressed psf file format is incorrect, please re-generate!\n");
4817 float psfVer = 0.0f;
4818 sscanf(buffer,
"FORMAT VERSION: %f\n", &psfVer);
4820 NAMD_die(
"The compressed psf file format is incorrect, please re-generate!\n");
4825 NAMD_die(
"UNABLE TO FIND NSEGMENTNAMES");
4826 sscanf(buffer,
"%d", &segNamePoolSize);
4828 if(segNamePoolSize!=0)
4830 for(
int i=0; i<segNamePoolSize; i++){
4832 sscanf(buffer,
"%s", strBuf);
4833 segNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4837 for(
int i=0; i<segNamePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4842 NAMD_die(
"UNABLE TO FIND NRESIDUENAMES");
4843 sscanf(buffer,
"%d", &resNamePoolSize);
4845 if(resNamePoolSize!=0)
4847 for(
int i=0; i<resNamePoolSize; i++){
4849 sscanf(buffer,
"%s", strBuf);
4850 resNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4854 for(
int i=0; i<resNamePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4859 NAMD_die(
"UNABLE TO FIND NATOMNAMES");
4860 sscanf(buffer,
"%d", &atomNamePoolSize);
4861 if(atomNamePoolSize!=0)
4863 for(
int i=0; i<atomNamePoolSize; i++){
4865 sscanf(buffer,
"%s", strBuf);
4866 atomNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4872 NAMD_die(
"UNABLE TO FIND NATOMTYPES");
4873 sscanf(buffer,
"%d", &atomTypePoolSize);
4875 if(atomTypePoolSize!=0)
4877 for(
int i=0; i<atomTypePoolSize; i++){
4879 sscanf(buffer,
"%s", strBuf);
4880 atomTypePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4884 for(
int i=0; i<atomTypePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4889 NAMD_die(
"UNABLE TO FIND NCHARGES");
4890 sscanf(buffer,
"%d", &chargePoolSize);
4891 if(chargePoolSize!=0)
4892 atomChargePool =
new Real[chargePoolSize];
4893 for(
int i=0; i<chargePoolSize; i++){
4895 sscanf(buffer,
"%f", atomChargePool+i);
4900 NAMD_die(
"UNABLE TO FIND NMASSES");
4901 sscanf(buffer,
"%d", &massPoolSize);
4903 atomMassPool =
new Real[massPoolSize];
4904 for(
int i=0; i<massPoolSize; i++){
4906 sscanf(buffer,
"%f", atomMassPool+i);
4911 NAMD_die(
"UNABLE TO FIND ATOMSIGS");
4912 sscanf(buffer,
"%d", &atomSigPoolSize);
4915 int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
4918 for(
int i=0; i<atomSigPoolSize; i++){
4922 NAMD_die(
"UNABLE TO FIND NBONDSIGS");
4923 sscanf(buffer,
"%d", &typeCnt);
4928 for(
int j=0; j<typeCnt; j++){
4930 sscanf(buffer,
"%d | %d | %d", &tmp1, &ttype, &tisReal);
4932 oneSig.offset[0] = tmp1;
4934 if(tisReal) numRealBonds++;
4940 NAMD_die(
"UNABLE TO FIND NTHETASIGS");
4941 sscanf(buffer,
"%d", &typeCnt);
4946 for(
int j=0; j<typeCnt; j++){
4948 sscanf(buffer,
"%d %d | %d | %d", &tmp1, &tmp2, &ttype, &tisReal);
4950 oneSig.offset[0] = tmp1;
4951 oneSig.offset[1] = tmp2;
4957 NAMD_die(
"UNABLE TO FIND NPHISIGS");
4958 sscanf(buffer,
"%d", &typeCnt);
4963 for(
int j=0; j<typeCnt; j++){
4965 sscanf(buffer,
"%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
4967 oneSig.offset[0] = tmp1;
4968 oneSig.offset[1] = tmp2;
4969 oneSig.offset[2] = tmp3;
4975 NAMD_die(
"UNABLE TO FIND NIMPHISIGS");
4976 sscanf(buffer,
"%d", &typeCnt);
4981 for(
int j=0; j<typeCnt; j++){
4983 sscanf(buffer,
"%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
4985 oneSig.offset[0] = tmp1;
4986 oneSig.offset[1] = tmp2;
4987 oneSig.offset[2] = tmp3;
4993 NAMD_die(
"UNABLE TO FIND NCRTERMSIGS");
4994 sscanf(buffer,
"%d", &typeCnt);
4999 for(
int j=0; j<typeCnt; j++){
5001 sscanf(buffer,
"%d %d %d %d %d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &ttype, &tisReal);
5003 oneSig.offset[0] = tmp1;
5004 oneSig.offset[1] = tmp2;
5005 oneSig.offset[2] = tmp3;
5006 oneSig.offset[3] = tmp4;
5007 oneSig.offset[4] = tmp5;
5008 oneSig.offset[5] = tmp6;
5009 oneSig.offset[6] = tmp7;
5017 NAMD_die(
"UNABLE TO FIND NEXCLSIGS");
5019 sscanf(buffer,
"%d", &exclSigPoolSize);
5021 vector<int> fullExcls;
5022 vector<int> modExcls;
5023 for(
int i=0; i<exclSigPoolSize; i++){
5025 for(
int j=0; j<fullExclCnt; j++)
5028 for(
int j=0; j<modExclCnt; j++)
5032 exclSigPool[i].setOffsets(fullExcls, modExcls);
5041 NAMD_die(
"UNABLE TO FIND NCLUSTERS");
5043 sscanf(buffer,
"%d", &numClusters);
5048 sscanf(buffer,
"%d", &numAtoms);
5052 NAMD_die(
"UNABLE TO FIND NHYDROGENGROUP");
5053 sscanf(buffer,
"%d", &numHydrogenGroups);
5057 NAMD_die(
"UNABLE TO FIND MAXHYDROGENGROUPSIZE");
5058 sscanf(buffer,
"%d", &maxHydrogenGroupSize);
5061 NAMD_die(
"UNABLE TO FIND NMIGRATIONGROUP");
5062 sscanf(buffer,
"%d", &numMigrationGroups);
5065 NAMD_die(
"UNABLE TO FIND MAXMIGRATIONGROUPSIZE");
5066 sscanf(buffer,
"%d", &maxMigrationGroupSize);
5068 int inputRigidType = -1;
5071 NAMD_die(
"UNABLE TO FIND RIGIDBONDTYPE");
5072 sscanf(buffer,
"%d", &inputRigidType);
5075 if(
simParams->rigidBonds != inputRigidType){
5076 const char *tmpstr[]={
"RIGID_NONE",
"RIGID_ALL",
"RIGID_WATER"};
5078 sprintf(errmsg,
"RIGIDBOND TYPE MISMATCH BETWEEN INPUT (%s) AND CURRENT RUN (%s)",
5079 tmpstr[inputRigidType], tmpstr[
simParams->rigidBonds]);
5087 NAMD_die(
"UNABLE TO FIND OCCUPANCYVALID");
5088 sscanf(buffer,
"%d", &isOccupancyValid);
5091 NAMD_die(
"UNABLE TO FIND TEMPFACTORVALID");
5092 sscanf(buffer,
"%d", &isBFactorValid);
5100 build_extra_bonds(params, cfgList->
find(
"extraBondsFile"));
5104 NAMD_die(
"UNABLE TO FIND DIHEDRALPARAMARRAY");
5113 NAMD_die(
"UNABLE TO FIND IMPROPERPARAMARRAY");
5130 void Molecule::read_binary_atom_info(
int fromAtomID,
int toAtomID,
InputAtomList& inAtoms){
5131 int numAtomsPar = toAtomID-fromAtomID+1;
5132 CmiAssert(numAtomsPar > 0);
5133 CmiAssert(inAtoms.
size() == numAtomsPar);
5152 eachAtomMass = NULL;
5153 eachAtomCharge = NULL;
5155 eachAtomExclSig = NULL;
5178 FILE *perAtomFile = fopen(
simParams->binAtomFile,
"rb");
5179 if (perAtomFile==NULL) {
5181 sprintf(err_msg,
"UNABLE TO OPEN THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s",
simParams->binAtomFile);
5187 flipNum((
char *)&rMagicNum,
sizeof(
int), 1);
5189 fread(&fMagicNum,
sizeof(
int), 1, perAtomFile);
5190 if (fMagicNum==magicNum) {
5192 }
else if (fMagicNum==rMagicNum) {
5196 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS CORRUPTED",
simParams->binAtomFile);
5200 float verNum = 0.0f;
5201 fread(&verNum,
sizeof(
float), 1, perAtomFile);
5202 if (needFlip)
flipNum((
char *)&verNum,
sizeof(
float), 1);
5205 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n",
simParams->binAtomFile);
5210 fread(&recSize,
sizeof(
int), 1, perAtomFile);
5211 if(needFlip)
flipNum((
char *)&recSize,
sizeof(
int), 1);
5214 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM RECORD SIZE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n",
simParams->binAtomFile);
5218 const int BUFELEMS = 32*1024;
5223 if ( _fseeki64(perAtomFile,startbyte,SEEK_CUR) )
5225 if ( fseeko(perAtomFile,startbyte,SEEK_CUR) )
5229 sprintf(errmsg,
"Error on seeking binary file %s",
simParams->binAtomFile);
5235 int atomsCnt = numAtomsPar;
5238 while(atomsCnt >= BUFELEMS) {
5239 if ( fread((
char *)elemsBuf,
sizeof(
OutputAtomRecord), BUFELEMS, perAtomFile) != BUFELEMS ) {
5241 sprintf(errmsg,
"Error on reading binary file %s",
simParams->binAtomFile);
5245 for(
int i=0; i<BUFELEMS; i++, curIdx++, oneRec++) {
5247 int aid = curIdx+fromAtomID;
5248 if(needFlip) oneRec->
flip();
5249 load_one_inputatom(aid, oneRec, fAtom);
5251 atomsCnt -= BUFELEMS;
5254 if ( fread(elemsBuf,
sizeof(
OutputAtomRecord), atomsCnt, perAtomFile) != atomsCnt ) {
5256 sprintf(errmsg,
"Error on reading binary file %s",
simParams->binAtomFile);
5260 for(
int i=curIdx; i<numAtomsPar; i++, oneRec++) {
5262 int aid = i+fromAtomID;
5263 if(needFlip) oneRec->
flip();
5264 load_one_inputatom(aid,oneRec,fAtom);
5267 if ( fclose(perAtomFile) ) {
5269 sprintf(errmsg,
"Error on closing binary file %s",
simParams->binAtomFile);
5278 is_atom_fixed(fromAtomID, &listIdx);
5279 for(
int i=listIdx; i<fixedAtomsSet->size(); i++){
5280 const AtomSet one = fixedAtomsSet->item(i);
5282 int sAtomId = one.aid1>fromAtomID ? one.aid1:fromAtomID;
5283 int eAtomId = one.aid2>toAtomID? toAtomID:one.aid2;
5284 for(
int j=sAtomId; j<=eAtomId; j++)
5285 inAtoms[j-fromAtomID].atomFixed = 1;
5292 char *thisAtomName = NULL;
5354 }
else if (thisAtomMass <= 0.05) {
5356 }
else if (thisAtomMass < 1.0) {
5358 }
else if (thisAtomMass <= 3.5) {
5360 }
else if (thisAtomName[0]==
'O' &&
5361 (thisAtomMass >= 14.0) && (thisAtomMass <= 18.0)) {
5375 void Molecule::build_excl_check_signatures(){
5377 for(
int i=0; i<exclSigPoolSize; i++){
5385 int fullMin, fullMax, modMin, modMax;
5393 if(fullMin < modMin)
5394 sigChk->
min = fullMin;
5396 sigChk->
min = modMin;
5397 if(fullMax < modMax)
5398 sigChk->
max = modMax;
5400 sigChk->
max = fullMax;
5408 iout <<
iWARN <<
"an empty exclusion signature with index " 5409 << i <<
"!\n" <<
endi;
5414 sigChk->
flags =
new char[sigChk->
max-sigChk->
min+1];
5415 memset(sigChk->
flags, 0,
sizeof(
char)*(sigChk->
max-sigChk->
min+1));
5436 void Molecule::load_atom_set(
StringList *setfile,
const char *setname,
5437 int *numAtomsInSet, AtomSetList **atomsSet)
const {
5438 if(setfile == NULL) {
5440 sprintf(errmsg,
"The text input file for %s atoms is not found!", setname);
5443 FILE *ifp = fopen(setfile->
data,
"r");
5447 sprintf(errmsg,
"ERROR IN OPENING %s ATOMS FILE: %s\n", setname, setfile->
data);
5452 int numLocalAtoms = 0;
5454 AtomSetList *localAtomsSet =
new AtomSetList();
5459 bool hasDash =
false;
5460 for(
int i=0; oneline[i] && i<128; i++){
5461 if(oneline[i]==
'-') {
5467 sscanf(oneline,
"%d-%d", &(one.aid1), &(one.aid2));
5468 if(one.aid1>one.aid2 || one.aid1<0 || one.aid2<0) {
5470 sprintf(errmsg,
"The input for %s atoms is wrong: %s\n", setname, oneline);
5473 numLocalAtoms += (one.aid2-one.aid1+1);
5475 sscanf(oneline,
"%d", &(one.aid1));
5478 sprintf(errmsg,
"The input for %s atoms is wrong: %s\n", setname, oneline);
5481 one.aid2 = one.aid1;
5484 localAtomsSet->add(one);
5488 std::sort(localAtomsSet->begin(), localAtomsSet->end());
5490 *numAtomsInSet = numLocalAtoms;
5491 *atomsSet = localAtomsSet;
5494 void Molecule::load_fixed_atoms(
StringList *fixedfile){
5495 load_atom_set(fixedfile,
"FIXED", &numFixedAtoms, &fixedAtomsSet);
5498 void Molecule::load_constrained_atoms(
StringList *constrainedfile){
5499 load_atom_set(constrainedfile,
"CONSTRAINED", &numConstraints, &constrainedAtomsSet);
5502 Bool Molecule::is_atom_in_set(AtomSetList *localAtomsSet,
int aid,
int *listIdx)
const {
5503 int idx = localAtomsSet->size();
5505 int lIdx = localAtomsSet->size()-1;
5507 while(rIdx <= lIdx){
5508 int mIdx = (rIdx+lIdx)/2;
5509 const AtomSet one = localAtomsSet->item(mIdx);
5515 }
else if(aid > one.aid1){
5519 if(listIdx) *listIdx = mIdx;
5526 if(listIdx) *listIdx = mIdx;
5532 if(listIdx) *listIdx = idx;
5550 #ifdef MEM_OPT_VERSION 5551 DebugM(3,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5560 <<
"******************************************\n" \
5561 <<
"NUM NAME TYPE RES MASS CHARGE CHARGE FEP-CHARGE" \
5562 <<
"SIGMA EPSILON SIGMA14 EPSILON14\n" \
5566 for (i=0; i<numAtoms; i++)
5568 const int vdw_type =
simParams->soluteScalingOn ?
5569 ((atoms[i].vdw_type >= LJtypecount) ?
5570 ss_vdw_type[atoms[i].vdw_type-LJtypecount] : atoms[i].vdw_type) : atoms[i].vdw_type;
5571 params->
get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, vdw_type);
5573 DebugM(3,i+1 <<
" " << atomNames[i].atomname \
5574 <<
" " << atomNames[i].atomtype <<
" " \
5575 << atomNames[i].resname <<
" " << atoms[i].mass \
5576 <<
" " << atoms[i].charge <<
" " << sigma \
5577 <<
" " << epsilon <<
" " << sigma14 \
5578 <<
" " << epsilon14 <<
"\n" \
5596 #ifdef MEM_OPT_VERSION 5597 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5603 DebugM(2,
"BOND LIST\n" <<
"********************************\n" \
5604 <<
"ATOM1 ATOM2 TYPE1 TYPE2 k x0" \
5607 for (i=0; i<numBonds; i++)
5611 DebugM(2,bonds[i].atom1+1 <<
" " \
5612 << bonds[i].atom2+1 <<
" " \
5613 << atomNames[bonds[i].atom1].atomtype <<
" " \
5614 << atomNames[bonds[i].atom2].atomtype <<
" " << k \
5615 <<
" " << x0 <<
endi);
5633 #ifdef MEM_OPT_VERSION 5634 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5638 DebugM(2,
"EXPLICIT EXCLUSION LIST\n" \
5639 <<
"********************************\n" \
5643 for (i=0; i<numExclusions; i++)
5645 DebugM(2,exclusions[i].atom1+1 <<
" " \
5646 << exclusions[i].atom2+1 <<
endi);
5663 #ifdef MEM_OPT_VERSION 5669 msg->
put(massPoolSize);
5670 msg->
put(massPoolSize, atomMassPool);
5672 msg->
put(chargePoolSize);
5673 msg->
put(chargePoolSize, atomChargePool);
5676 msg->
put(atomSigPoolSize);
5677 for(
int i=0; i<atomSigPoolSize; i++)
5681 msg->
put(exclSigPoolSize);
5682 for(
int i=0; i<exclSigPoolSize; i++)
5683 exclSigPool[i].pack(msg);
5685 msg->
put(numHydrogenGroups);
5686 msg->
put(maxHydrogenGroupSize);
5687 msg->
put(numMigrationGroups);
5688 msg->
put(maxMigrationGroupSize);
5689 msg->
put(isOccupancyValid);
5690 msg->
put(isBFactorValid);
5693 msg->
put(atomNamePoolSize);
5694 for(
int i=0; i<atomNamePoolSize;i++) {
5701 int numFixedAtomsSet = fixedAtomsSet->size();
5702 msg->
put(numFixedAtoms);
5703 msg->
put(numFixedAtomsSet);
5704 msg->
put(numFixedAtomsSet*
sizeof(AtomSet), (
char *)(fixedAtomsSet->begin()));
5708 int numConstrainedAtomsSet = constrainedAtomsSet->size();
5709 msg->
put(numConstraints);
5710 msg->
put(numConstrainedAtomsSet);
5711 msg->
put(numConstrainedAtomsSet*
sizeof(AtomSet), (
char *)(constrainedAtomsSet->begin()));
5716 msg->
put(numAtoms*
sizeof(
Atom), (
char*)atoms);
5719 msg->
put(numRealBonds);
5724 msg->
put(numBonds*
sizeof(
Bond), (
char*)bonds);
5728 msg->
put(numAngles);
5731 msg->
put(numAngles*
sizeof(
Angle), (
char*)angles);
5735 msg->
put(numDihedrals);
5738 msg->
put(numDihedrals*
sizeof(
Dihedral), (
char*)dihedrals);
5742 msg->
put(num_alch_unpert_Bonds);
5743 msg->
put(num_alch_unpert_Bonds*
sizeof(
Bond), (
char*)alch_unpert_bonds);
5745 msg->
put(num_alch_unpert_Angles);
5746 msg->
put(num_alch_unpert_Angles*
sizeof(
Angle), (
char*)alch_unpert_angles);
5748 msg->
put(num_alch_unpert_Dihedrals);
5749 msg->
put(num_alch_unpert_Dihedrals*
sizeof(
Dihedral), (
char*)alch_unpert_dihedrals);
5753 msg->
put(numImpropers);
5756 msg->
put(numImpropers*
sizeof(
Improper), (
char*)impropers);
5760 msg->
put(numCrossterms);
5763 msg->
put(numCrossterms*
sizeof(
Crossterm), (
char*)crossterms);
5767 msg->
put(numDonors);
5770 msg->
put(numDonors*
sizeof(
Bond), (
char*)donors);
5774 msg->
put(numAcceptors);
5777 msg->
put(numAcceptors*
sizeof(
Bond), (
char*)acceptors);
5781 msg->
put(numExclusions);
5784 msg->
put(numExclusions*
sizeof(
Exclusion), (
char*)exclusions);
5789 msg->
put(numConstraints);
5791 msg->
put(numAtoms, consIndexes);
5795 msg->
put(numConstraints*
sizeof(ConstraintParams), (
char*)consParams);
5804 DebugM(3,
"Sending gridforce info\n" <<
endi);
5805 msg->
put(numGridforceGrids);
5807 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
5808 msg->
put(numGridforces[gridnum]);
5809 msg->
put(numAtoms, gridfrcIndexes[gridnum]);
5810 if (numGridforces[gridnum])
5812 msg->
put(numGridforces[gridnum]*
sizeof(GridforceParams), (
char*)gridfrcParams[gridnum]);
5823 msg->
put(numStirredAtoms);
5825 msg->
put(numAtoms, stirIndexes);
5827 if (numStirredAtoms)
5830 msg->
put(numStirredAtoms*
sizeof(StirParams), (
char*)stirParams);
5837 msg->
put(numMovDrag);
5838 msg->
put(numAtoms, movDragIndexes);
5841 msg->
put(numMovDrag*
sizeof(MovDragParams), (
char*)movDragParams);
5847 msg->
put(numRotDrag);
5848 msg->
put(numAtoms, rotDragIndexes);
5851 msg->
put(numRotDrag*
sizeof(RotDragParams), (
char*)rotDragParams);
5857 msg->
put(numConsTorque);
5858 msg->
put(numAtoms, consTorqueIndexes);
5861 msg->
put(numConsTorque*
sizeof(ConsTorqueParams), (
char*)consTorqueParams);
5867 { msg->
put(numConsForce);
5868 msg->
put(numAtoms, consForceIndexes);
5870 msg->
put(numConsForce*
sizeof(
Vector), (
char*)consForce);
5875 msg->
put(numMolecules);
5876 msg->
put(numLargeMolecules);
5877 msg->
put(numAtoms, moleculeAtom);
5878 msg->
put(numMolecules+1, moleculeStartIndex);
5882 msg->
put(numExPressureAtoms);
5883 msg->
put(numAtoms, exPressureAtomFlags);
5886 #ifndef MEM_OPT_VERSION 5890 msg->
put(numAtoms, langevinParams);
5896 msg->
put(numFixedAtoms);
5897 msg->
put(numAtoms, fixedAtomFlags);
5898 msg->
put(numFixedRigidBonds);
5903 msg->
put(numAtoms, qmAtomGroup);
5904 msg->
put(qmNumQMAtoms);
5905 msg->
put(qmNumQMAtoms, qmAtmChrg);
5906 msg->
put(qmNumQMAtoms, qmAtmIndx);
5908 msg->
put(qmNumBonds);
5909 msg->
put(qmMeNumBonds);
5910 msg->
put(qmMeNumBonds, qmMeMMindx);
5911 msg->
put(qmMeNumBonds, qmMeQMGrp);
5913 msg->
put(qmNumGrps);
5914 msg->
put(qmNumGrps, qmGrpID);
5915 msg->
put(qmNumGrps, qmCustPCSizes);
5916 msg->
put(qmTotCustPCs);
5917 msg->
put(qmTotCustPCs, qmCustomPCIdxs);
5923 msg->
put(numFepInitial);
5924 msg->
put(numFepFinal);
5925 msg->
put(numAtoms*
sizeof(
char), (
char*)fepAtomFlags);
5930 msg->
put(numAtoms*
sizeof(
char), (
char*)ssAtomFlags);
5931 msg->
put(ss_num_vdw_params);
5933 msg->
put(numAtoms*
sizeof(
int), (
char*)ss_index);
5936 #ifdef OPENATOM_VERSION 5939 msg->
put(numFepInitial);
5940 msg->
put(numAtoms*
sizeof(
char), (
char*)fepAtomFlags);
5942 #endif //OPENATOM_VERSION 5945 msg->
put(is_lonepairs_psf);
5946 if (is_lonepairs_psf) {
5947 msg->
put(numLphosts);
5948 msg->
put(numLphosts*
sizeof(
Lphost), (
char*)lphosts);
5950 msg->
put(is_drude_psf);
5953 msg->
put(numTholes);
5954 msg->
put(numTholes*
sizeof(
Thole), (
char*)tholes);
5955 msg->
put(numAnisos);
5956 msg->
put(numAnisos*
sizeof(
Aniso), (
char*)anisos);
5958 msg->
put(numZeroMassAtoms);
5963 msg->
put(numAtoms, (
int*)lcpoParamType);
5974 msg->
put((numAtoms),pointerToLJBeg);
5975 msg->
put((numAtoms),pointerToLJEnd);
5985 msg->
put((numAtoms),pointerToGaussBeg);
5986 msg->
put((numAtoms),pointerToGaussEnd);
5994 #ifdef MEM_OPT_VERSION 5996 build_excl_check_signatures();
5999 numBonds = numCalcBonds = 0;
6000 numAngles = numCalcAngles = 0;
6001 numDihedrals = numCalcDihedrals = 0;
6002 numImpropers = numCalcImpropers = 0;
6003 numCrossterms = numCalcCrossterms = 0;
6004 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;
6012 build_lists_by_atom();
6032 #ifdef MEM_OPT_VERSION 6036 msg->
get(massPoolSize);
6037 if(atomMassPool)
delete [] atomMassPool;
6038 atomMassPool =
new Real[massPoolSize];
6039 msg->
get(massPoolSize, atomMassPool);
6041 msg->
get(chargePoolSize);
6042 if(atomChargePool)
delete [] atomChargePool;
6043 atomChargePool =
new Real[chargePoolSize];
6044 msg->
get(chargePoolSize, atomChargePool);
6047 msg->
get(atomSigPoolSize);
6050 for(
int i=0; i<atomSigPoolSize; i++)
6054 msg->
get(exclSigPoolSize);
6055 if(exclSigPool)
delete [] exclSigPool;
6057 for(
int i=0; i<exclSigPoolSize; i++)
6058 exclSigPool[i].unpack(msg);
6060 msg->
get(numHydrogenGroups);
6061 msg->
get(maxHydrogenGroupSize);
6062 msg->
get(numMigrationGroups);
6063 msg->
get(maxMigrationGroupSize);
6064 msg->
get(isOccupancyValid);
6065 msg->
get(isBFactorValid);
6068 msg->
get(atomNamePoolSize);
6070 for(
int i=0; i<atomNamePoolSize;i++) {
6078 int numFixedAtomsSet;
6079 msg->
get(numFixedAtoms);
6080 msg->
get(numFixedAtomsSet);
6081 fixedAtomsSet =
new AtomSetList(numFixedAtomsSet);
6082 msg->
get(numFixedAtomsSet*
sizeof(AtomSet), (
char *)(fixedAtomsSet->begin()));
6086 int numConstrainedAtomsSet;
6087 msg->
get(numConstraints);
6088 msg->
get(numConstrainedAtomsSet);
6089 constrainedAtomsSet =
new AtomSetList(numConstrainedAtomsSet);
6090 msg->
get(numConstrainedAtomsSet*
sizeof(AtomSet), (
char *)(constrainedAtomsSet->begin()));
6095 atoms=
new Atom[numAtoms];
6096 msg->
get(numAtoms*
sizeof(
Atom), (
char*)atoms);
6099 msg->
get(numRealBonds);
6104 bonds=
new Bond[numBonds];
6105 msg->
get(numBonds*
sizeof(
Bond), (
char*)bonds);
6109 msg->
get(numAngles);
6113 angles=
new Angle[numAngles];
6114 msg->
get(numAngles*
sizeof(
Angle), (
char*)angles);
6118 msg->
get(numDihedrals);
6121 delete [] dihedrals;
6122 dihedrals=
new Dihedral[numDihedrals];
6123 msg->
get(numDihedrals*
sizeof(
Dihedral), (
char*)dihedrals);
6127 msg->
get(num_alch_unpert_Bonds);
6128 alch_unpert_bonds=
new Bond[num_alch_unpert_Bonds];
6129 msg->
get(num_alch_unpert_Bonds*
sizeof(
Bond), (
char*)alch_unpert_bonds);
6131 msg->
get(num_alch_unpert_Angles);
6132 alch_unpert_angles=
new Angle[num_alch_unpert_Angles];
6133 msg->
get(num_alch_unpert_Angles*
sizeof(
Angle), (
char*)alch_unpert_angles);
6135 msg->
get(num_alch_unpert_Dihedrals);
6136 alch_unpert_dihedrals=
new Dihedral[num_alch_unpert_Dihedrals];
6137 msg->
get(num_alch_unpert_Dihedrals*
sizeof(
Dihedral), (
char*)alch_unpert_dihedrals);
6141 msg->
get(numImpropers);
6144 delete [] impropers;
6145 impropers=
new Improper[numImpropers];
6146 msg->
get(numImpropers*
sizeof(
Improper), (
char*)impropers);
6150 msg->
get(numCrossterms);
6153 delete [] crossterms;
6154 crossterms=
new Crossterm[numCrossterms];
6155 msg->
get(numCrossterms*
sizeof(
Crossterm), (
char*)crossterms);
6159 msg->
get(numDonors);
6163 donors=
new Bond[numDonors];
6164 msg->
get(numDonors*
sizeof(
Bond), (
char*)donors);
6168 msg->
get(numAcceptors);
6171 delete [] acceptors;
6172 acceptors=
new Bond[numAcceptors];
6173 msg->
get(numAcceptors*
sizeof(
Bond), (
char*)acceptors);
6177 msg->
get(numExclusions);
6180 delete [] exclusions;
6181 exclusions=
new Exclusion[numExclusions];
6182 msg->
get(numExclusions*
sizeof(
Exclusion), (
char*)exclusions);
6188 msg->
get(numConstraints);
6190 delete [] consIndexes;
6191 consIndexes =
new int32[numAtoms];
6193 msg->
get(numAtoms, consIndexes);
6197 delete [] consParams;
6198 consParams =
new ConstraintParams[numConstraints];
6200 msg->
get(numConstraints*
sizeof(ConstraintParams), (
char*)consParams);
6208 DebugM(3,
"Receiving gridforce info\n");
6210 msg->
get(numGridforceGrids);
6212 DebugM(3,
"numGridforceGrids = " << numGridforceGrids <<
"\n");
6214 delete [] numGridforces;
6215 numGridforces =
new int[numGridforceGrids];
6217 delete [] gridfrcIndexes;
6218 delete [] gridfrcParams;
6219 delete [] gridfrcGrid;
6220 gridfrcIndexes =
new int32*[numGridforceGrids];
6221 gridfrcParams =
new GridforceParams*[numGridforceGrids];
6224 int grandTotalGrids = 0;
6225 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6226 msg->
get(numGridforces[gridnum]);
6228 gridfrcIndexes[gridnum] =
new int32[numAtoms];
6229 msg->
get(numAtoms, gridfrcIndexes[gridnum]);
6231 if (numGridforces[gridnum])
6233 gridfrcParams[gridnum] =
new GridforceParams[numGridforces[gridnum]];
6234 msg->
get(numGridforces[gridnum]*
sizeof(GridforceParams), (
char*)gridfrcParams[gridnum]);
6247 msg->
get(numStirredAtoms);
6249 delete [] stirIndexes;
6250 stirIndexes =
new int32[numAtoms];
6252 msg->
get(numAtoms, stirIndexes);
6254 if (numStirredAtoms)
6256 delete [] stirParams;
6257 stirParams =
new StirParams[numStirredAtoms];
6259 msg->
get(numStirredAtoms*
sizeof(StirParams), (
char*)stirParams);
6265 msg->
get(numMovDrag);
6266 delete [] movDragIndexes;
6267 movDragIndexes =
new int32[numAtoms];
6268 msg->
get(numAtoms, movDragIndexes);
6271 delete [] movDragParams;
6272 movDragParams =
new MovDragParams[numMovDrag];
6273 msg->
get(numMovDrag*
sizeof(MovDragParams), (
char*)movDragParams);
6279 msg->
get(numRotDrag);
6280 delete [] rotDragIndexes;
6281 rotDragIndexes =
new int32[numAtoms];
6282 msg->
get(numAtoms, rotDragIndexes);
6285 delete [] rotDragParams;
6286 rotDragParams =
new RotDragParams[numRotDrag];
6287 msg->
get(numRotDrag*
sizeof(RotDragParams), (
char*)rotDragParams);
6293 msg->
get(numConsTorque);
6294 delete [] consTorqueIndexes;
6295 consTorqueIndexes =
new int32[numAtoms];
6296 msg->
get(numAtoms, consTorqueIndexes);
6299 delete [] consTorqueParams;
6300 consTorqueParams =
new ConsTorqueParams[numConsTorque];
6301 msg->
get(numConsTorque*
sizeof(ConsTorqueParams), (
char*)consTorqueParams);
6307 { msg->
get(numConsForce);
6308 delete [] consForceIndexes;
6309 consForceIndexes =
new int32[numAtoms];
6310 msg->
get(numAtoms, consForceIndexes);
6312 {
delete [] consForce;
6313 consForce =
new Vector[numConsForce];
6314 msg->
get(numConsForce*
sizeof(
Vector), (
char*)consForce);
6320 msg->
get(numMolecules);
6321 msg->
get(numLargeMolecules);
6322 delete [] moleculeAtom;
6323 delete [] moleculeStartIndex;
6324 moleculeAtom =
new int32[numAtoms];
6325 moleculeStartIndex =
new int32[numMolecules+1];
6326 msg->
get(numAtoms, moleculeAtom);
6327 msg->
get(numMolecules+1, moleculeStartIndex);
6331 exPressureAtomFlags =
new int32[numAtoms];
6332 msg->
get(numExPressureAtoms);
6333 msg->
get(numAtoms, exPressureAtomFlags);
6336 #ifndef MEM_OPT_VERSION 6340 delete [] langevinParams;
6341 langevinParams =
new Real[numAtoms];
6343 msg->
get(numAtoms, langevinParams);
6349 delete [] fixedAtomFlags;
6350 fixedAtomFlags =
new int32[numAtoms];
6352 msg->
get(numFixedAtoms);
6353 msg->
get(numAtoms, fixedAtomFlags);
6354 msg->
get(numFixedRigidBonds);
6359 if( qmAtomGroup != 0)
6360 delete [] qmAtomGroup;
6361 qmAtomGroup =
new Real[numAtoms];
6363 msg->
get(numAtoms, qmAtomGroup);
6365 msg->
get(qmNumQMAtoms);
6368 delete [] qmAtmChrg;
6369 qmAtmChrg =
new Real[qmNumQMAtoms];
6371 msg->
get(qmNumQMAtoms, qmAtmChrg);
6374 delete [] qmAtmIndx;
6375 qmAtmIndx =
new int[qmNumQMAtoms];
6377 msg->
get(qmNumQMAtoms, qmAtmIndx);
6381 msg->
get(qmNumBonds);
6383 msg->
get(qmMeNumBonds);
6385 if( qmMeMMindx != 0)
6386 delete [] qmMeMMindx;
6387 qmMeMMindx =
new int[qmMeNumBonds];
6389 msg->
get(qmMeNumBonds, qmMeMMindx);
6392 delete [] qmMeQMGrp;
6393 qmMeQMGrp =
new Real[qmMeNumBonds];
6395 msg->
get(qmMeNumBonds, qmMeQMGrp);
6399 msg->
get(qmNumGrps);
6403 qmGrpID =
new Real[qmNumGrps];
6404 msg->
get(qmNumGrps, qmGrpID);
6406 if( qmCustPCSizes != 0)
6407 delete [] qmCustPCSizes;
6408 qmCustPCSizes =
new int[qmNumGrps];
6409 msg->
get(qmNumGrps, qmCustPCSizes);
6411 msg->
get(qmTotCustPCs);
6413 if( qmCustomPCIdxs != 0)
6414 delete [] qmCustomPCIdxs;
6415 qmCustomPCIdxs =
new int[qmTotCustPCs];
6416 msg->
get(qmTotCustPCs, qmCustomPCIdxs);
6422 delete [] fepAtomFlags;
6423 fepAtomFlags =
new unsigned char[numAtoms];
6425 msg->
get(numFepInitial);
6426 msg->
get(numFepFinal);
6427 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)fepAtomFlags);
6433 delete [] ssAtomFlags;
6434 delete [] ss_vdw_type;
6436 ssAtomFlags =
new unsigned char[numAtoms];
6438 ss_index =
new int [numAtoms];
6439 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)ssAtomFlags);
6440 msg->
get(ss_num_vdw_params);
6442 msg->
get(numAtoms*
sizeof(
int), (
char*)ss_index);
6445 #ifdef OPENATOM_VERSION 6448 delete [] fepAtomFlags;
6449 fepAtomFlags =
new unsigned char[numAtoms];
6451 msg->
get(numFepInitial);
6452 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)fepAtomFlags);
6453 #endif //OPENATOM_VERSION 6456 msg->
get(is_lonepairs_psf);
6457 if (is_lonepairs_psf) {
6458 msg->
get(numLphosts);
6460 lphosts =
new Lphost[numLphosts];
6461 msg->
get(numLphosts*
sizeof(
Lphost), (
char*)lphosts);
6463 msg->
get(is_drude_psf);
6465 delete[] drudeConsts;
6468 msg->
get(numTholes);
6470 tholes =
new Thole[numTholes];
6471 msg->
get(numTholes*
sizeof(
Thole), (
char*)tholes);
6472 msg->
get(numAnisos);
6474 anisos =
new Aniso[numAnisos];
6475 msg->
get(numAnisos*
sizeof(
Aniso), (
char*)anisos);
6477 msg->
get(numZeroMassAtoms);
6482 delete [] lcpoParamType;
6483 lcpoParamType =
new int[numAtoms];
6484 msg->
get(numAtoms, (
int*)lcpoParamType);
6503 delete [] gromacsPair_type;
6506 delete [] pointerToLJBeg;
6507 pointerToLJBeg =
new int[numAtoms];
6508 msg->
get((numAtoms),pointerToLJBeg);
6509 delete [] pointerToLJEnd;
6510 pointerToLJEnd =
new int[numAtoms];
6511 msg->
get((numAtoms),pointerToLJEnd);
6524 delete [] indxGaussA;
6527 delete [] indxGaussB;
6545 delete [] gRepulsive;
6548 delete [] pointerToGaussBeg;
6549 pointerToGaussBeg =
new int[numAtoms];
6550 msg->
get((numAtoms),pointerToGaussBeg);
6551 delete [] pointerToGaussEnd;
6552 pointerToGaussEnd =
new int[numAtoms];
6553 msg->
get((numAtoms),pointerToGaussEnd);
6561 #ifdef MEM_OPT_VERSION 6563 build_excl_check_signatures();
6566 numBonds = numCalcBonds = 0;
6567 numAngles = numCalcAngles = 0;
6568 numDihedrals = numCalcDihedrals = 0;
6569 numImpropers = numCalcImpropers = 0;
6570 numCrossterms = numCalcCrossterms = 0;
6571 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;
6579 build_atom_status();
6580 build_lists_by_atom();
6619 DebugM(3,
"Entered build_gridforce_params multi...\n");
6624 numGridforceGrids = 0;
6625 while (mgridParams != NULL) {
6626 numGridforceGrids++;
6627 mgridParams = mgridParams->
next;
6630 DebugM(3,
"numGridforceGrids = " << numGridforceGrids <<
"\n");
6631 gridfrcIndexes =
new int32*[numGridforceGrids];
6632 gridfrcParams =
new GridforceParams*[numGridforceGrids];
6634 numGridforces =
new int[numGridforceGrids];
6636 int grandTotalGrids = 0;
6638 mgridParams =
simParams->mgridforcelist.get_first();
6639 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6640 int current_index=0;
6647 if (mgridParams == NULL) {
6648 NAMD_die(
"Problem with mgridParams!");
6654 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, gridforceFile required.");
6661 if ( (cwd == NULL) || (mgridParams->
gridforceFile[0] ==
'/') )
6667 strcpy(filename, cwd);
6671 kPDB =
new PDB(filename);
6674 NAMD_die(
"Memory allocation failed in Molecule::build_gridforce_params");
6679 NAMD_die(
"Number of atoms in grid force PDB doesn't match coordinate PDB");
6698 else if (strcasecmp(mgridParams->
gridforceCol,
"Y") == 0)
6702 else if (strcasecmp(mgridParams->
gridforceCol,
"Z") == 0)
6706 else if (strcasecmp(mgridParams->
gridforceCol,
"O") == 0)
6710 else if (strcasecmp(mgridParams->
gridforceCol,
"B") == 0)
6716 NAMD_die(
"gridforcecol must have value of X, Y, Z, O, or B");
6749 NAMD_die(
"gridforcechargecol must have value of X, Y, Z, O, or B");
6754 NAMD_die(
"gridforcecol and gridforcechargecol cannot have same value");
6761 gridfrcIndexes[gridnum] =
new int32[numAtoms];
6763 if (gridfrcIndexes[gridnum] == NULL)
6765 NAMD_die(
"memory allocation failed in Molecule::build_gridforce_params()");
6769 for (i=0; i<numAtoms; i++)
6775 kval = (kPDB->
atom(i))->xcoor();
6778 kval = (kPDB->
atom(i))->ycoor();
6781 kval = (kPDB->
atom(i))->zcoor();
6784 kval = (kPDB->
atom(i))->occupancy();
6787 kval = (kPDB->
atom(i))->temperaturefactor();
6794 gridfrcIndexes[gridnum][i] = current_index;
6800 gridfrcIndexes[gridnum][i] = -1;
6804 if (current_index == 0)
6807 iout <<
iWARN <<
"NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" <<
endi;
6812 gridfrcParams[gridnum] =
new GridforceParams[current_index];
6813 if (gridfrcParams[gridnum] == NULL)
6815 NAMD_die(
"memory allocation failed in Molecule::build_gridforce_params");
6819 numGridforces[gridnum] = current_index;
6823 for (i=0; i<numAtoms; i++)
6825 if (gridfrcIndexes[gridnum][i] != -1)
6831 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->xcoor();
6834 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->ycoor();
6837 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->zcoor();
6840 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->occupancy();
6843 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->temperaturefactor();
6851 #ifdef MEM_OPT_VERSION 6852 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
6854 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
6858 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->xcoor();
6861 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->ycoor();
6864 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->zcoor();
6867 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->occupancy();
6870 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->temperaturefactor();
6891 strcpy(potfilename, cwd);
6898 DebugM(3,
"allocating GridforceGrid(" << gridnum <<
")\n");
6902 DebugM(4,
"grandTotalGrids = " << grandTotalGrids <<
"\n" <<
endi);
6905 mgridParams = mgridParams->
next;
6910 #ifdef MEM_OPT_VERSION 6911 void Molecule::delEachAtomSigs(){
6917 if(CmiMyRank())
return;
6919 delete [] eachAtomSig;
6920 delete [] eachAtomExclSig;
6922 eachAtomExclSig = NULL;
6925 void Molecule::delChargeSpace(){
6926 if(CmiMyRank())
return;
6928 delete [] atomChargePool;
6929 delete [] eachAtomCharge;
6930 atomChargePool = NULL;
6931 eachAtomCharge = NULL;
6934 void Molecule::delMassSpace(){
6935 if(CmiMyRank())
return;
6937 delete [] atomMassPool;
6938 delete [] eachAtomMass;
6939 atomMassPool = NULL;
6940 eachAtomMass = NULL;
6943 void Molecule::delClusterSigs() {
6944 if(CmiMyRank())
return;
6946 delete [] clusterSigs;
6950 void Molecule::delAtomNames(){
6951 if(CmiMyRank())
return;
6953 delete [] atomNames;
6958 void Molecule::delFixedAtoms(){
6959 if(CmiMyRank())
return;
6960 delete fixedAtomsSet;
6961 fixedAtomsSet = NULL;
6966 #endif // MOLECULE2_C undefined = first object file 6967 #ifdef MOLECULE2_C // second object file 6999 int current_index=0;
7009 if (consref == NULL)
7011 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, consref required.");
7012 refPDB = initial_pdb;
7016 if (consref->
next != NULL)
7018 NAMD_die(
"Multiple definitions of constraint reference file in configruation file");
7021 if ( (cwd == NULL) || (consref->
data[0] ==
'/') )
7023 strcpy(filename, consref->
data);
7027 strcpy(filename, cwd);
7028 strcat(filename, consref->
data);
7031 refPDB =
new PDB(filename);
7032 if ( refPDB == NULL )
7034 NAMD_die(
"Memory allocation failed in Molecule::build_constraint_params");
7039 NAMD_die(
"Number of atoms in constraint reference PDB doesn't match coordinate PDB");
7046 if (conskfile == NULL)
7048 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, conskfile required.");
7053 if (conskfile->
next != NULL)
7055 NAMD_die(
"Multiple definitions of constraint constant file in configuration file");
7058 if ( (consref != NULL) && (strcasecmp(consref->
data, conskfile->
data) == 0) )
7065 if ( (cwd == NULL) || (conskfile->
data[0] ==
'/') )
7067 strcpy(filename, conskfile->
data);
7071 strcpy(filename, cwd);
7072 strcat(filename, conskfile->
data);
7075 kPDB =
new PDB(filename);
7078 NAMD_die(
"Memory allocation failed in Molecule::build_constraint_params");
7083 NAMD_die(
"Number of atoms in constraint constant PDB doesn't match coordinate PDB");
7093 if (conskcol == NULL)
7099 if (conskcol->
next != NULL)
7101 NAMD_die(
"Multiple definitions of harmonic constraint column in config file");
7104 if (strcasecmp(conskcol->
data,
"X") == 0)
7108 else if (strcasecmp(conskcol->
data,
"Y") == 0)
7112 else if (strcasecmp(conskcol->
data,
"Z") == 0)
7116 else if (strcasecmp(conskcol->
data,
"O") == 0)
7120 else if (strcasecmp(conskcol->
data,
"B") == 0)
7126 NAMD_die(
"conskcol must have value of X, Y, Z, O, or B");
7133 consIndexes =
new int32[numAtoms];
7135 if (consIndexes == NULL)
7137 NAMD_die(
"memory allocation failed in Molecule::build_constraint_params()");
7141 for (i=0; i<numAtoms; i++)
7147 kval = (kPDB->
atom(i))->xcoor();
7150 kval = (kPDB->
atom(i))->ycoor();
7153 kval = (kPDB->
atom(i))->zcoor();
7156 kval = (kPDB->
atom(i))->occupancy();
7159 kval = (kPDB->
atom(i))->temperaturefactor();
7166 consIndexes[i] = current_index;
7172 consIndexes[i] = -1;
7176 if (current_index == 0)
7179 iout <<
iWARN <<
"NO CONSTRAINED ATOMS WERE FOUND, BUT CONSTRAINTS ARE ON . . .\n" <<
endi;
7184 consParams =
new ConstraintParams[current_index];
7186 if (consParams == NULL)
7188 NAMD_die(
"memory allocation failed in Molecule::build_constraint_params");
7192 numConstraints = current_index;
7196 for (i=0; i<numAtoms; i++)
7198 if (consIndexes[i] != -1)
7204 consParams[consIndexes[i]].k = (kPDB->
atom(i))->xcoor();
7207 consParams[consIndexes[i]].k = (kPDB->
atom(i))->ycoor();
7210 consParams[consIndexes[i]].k = (kPDB->
atom(i))->zcoor();
7213 consParams[consIndexes[i]].k = (kPDB->
atom(i))->occupancy();
7216 consParams[consIndexes[i]].k = (kPDB->
atom(i))->temperaturefactor();
7221 consParams[consIndexes[i]].refPos.x = (refPDB->
atom(i))->xcoor();
7222 consParams[consIndexes[i]].refPos.y = (refPDB->
atom(i))->ycoor();
7223 consParams[consIndexes[i]].refPos.z = (refPDB->
atom(i))->zcoor();
7228 if (consref != NULL)
7233 if ((conskfile != NULL) &&
7234 !((consref != NULL) &&
7235 (strcasecmp(consref->
data, conskfile->
data) == 0)
7274 int current_index=0;
7283 if (movDragFile == NULL) {
7284 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, movDragFile required.");
7289 if (movDragFile->
next != NULL) {
7290 NAMD_die(
"Multiple definitions of moving drag tag file in configuration file");
7293 if ( (cwd == NULL) || (movDragFile->
data[0] ==
'/') ) {
7294 strcpy(mainfilename, movDragFile->
data);
7296 strcpy(mainfilename, cwd);
7297 strcat(mainfilename, movDragFile->
data);
7300 tPDB =
new PDB(mainfilename);
7301 if ( tPDB == NULL ) {
7302 NAMD_die(
"Memory allocation failed in Molecule::build_movdrag_params");
7306 NAMD_die(
"Number of atoms in moving drag tag PDB doesn't match coordinate PDB");
7314 if (movDragVelFile == NULL) {
7315 if (movDragFile == NULL) {
7316 NAMD_die(
"Moving drag velocity file can not be same as coordinate PDB file");
7318 if (movDragVelFile->
next != NULL) {
7319 NAMD_die(
"Multiple definitions of moving drag velocity file in configuration file");
7326 if ( (cwd == NULL) || (movDragVelFile->
data[0] ==
'/') ) {
7327 strcpy(velfilename, movDragVelFile->
data);
7329 strcpy(velfilename, cwd);
7330 strcat(velfilename, movDragVelFile->
data);
7333 vPDB =
new PDB(velfilename);
7334 if ( vPDB == NULL ) {
7335 NAMD_die(
"Memory allocation failed in Molecule::build_movdrag_params");
7339 NAMD_die(
"Number of atoms in moving drag velocity PDB doesn't match coordinate PDB");
7351 if (movDragCol == NULL) {
7354 if (movDragCol->
next != NULL) {
7355 NAMD_die(
"Multiple definitions of drag column in config file");
7358 if (movDragFile == NULL
7359 && strcasecmp(movDragCol->
data,
"B")
7360 && strcasecmp(movDragCol->
data,
"O")) {
7361 NAMD_die(
"Can not read moving drag tags from X, Y, or Z column of the coordinate or velocity file");
7363 if (!strcasecmp(movDragCol->
data,
"X")) {
7365 }
else if (!strcasecmp(movDragCol->
data,
"Y")) {
7367 }
else if (!strcasecmp(movDragCol->
data,
"Z")) {
7369 }
else if (!strcasecmp(movDragCol->
data,
"O")) {
7371 }
else if (!strcasecmp(movDragCol->
data,
"B")) {
7375 NAMD_die(
"movDragCol must have value of X, Y, Z, O, or B");
7382 movDragIndexes =
new int32[numAtoms];
7383 if (movDragIndexes == NULL) {
7384 NAMD_die(
"memory allocation failed in Molecule::build_movdrag_params()");
7388 for (i=0; i<numAtoms; i++) {
7391 dtval = (tPDB->
atom(i))->xcoor();
7394 dtval = (tPDB->
atom(i))->ycoor();
7397 dtval = (tPDB->
atom(i))->zcoor();
7400 dtval = (tPDB->
atom(i))->occupancy();
7403 dtval = (tPDB->
atom(i))->temperaturefactor();
7409 movDragIndexes[i] = current_index;
7413 movDragIndexes[i] = -1;
7417 if (current_index == 0) {
7419 iout <<
iWARN <<
"NO DRAGGED ATOMS WERE FOUND, BUT MOVING DRAG IS ON . . . " <<
endi;
7422 movDragParams =
new MovDragParams[current_index];
7423 if (movDragParams == NULL) {
7424 NAMD_die(
"memory allocation failed in Molecule::build_movdrag_params");
7428 numMovDrag = current_index;
7432 for (i=0; i<numAtoms; i++) {
7433 if (movDragIndexes[i] != -1) {
7434 movDragParams[movDragIndexes[i]].v[0] = (vPDB->
atom(i))->xcoor();
7435 movDragParams[movDragIndexes[i]].v[1] = (vPDB->
atom(i))->ycoor();
7436 movDragParams[movDragIndexes[i]].v[2] = (vPDB->
atom(i))->zcoor();
7440 if (movDragFile != NULL)
delete tPDB;
7441 if (movDragVelFile != NULL)
delete vPDB;
7478 PDB *tPDB, *aPDB, *pPDB, *vPDB;
7480 int current_index=0;
7493 if (rotDragFile == NULL) {
7494 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, rotDragFile required.");
7499 if (rotDragFile->
next != NULL) {
7500 NAMD_die(
"Multiple definitions of rotating drag tag file in configuration file");
7503 if ( (cwd == NULL) || (rotDragFile->
data[0] ==
'/') ) {
7504 strcpy(mainfilename, rotDragFile->
data);
7506 strcpy(mainfilename, cwd);
7507 strcat(mainfilename, rotDragFile->
data);
7510 tPDB =
new PDB(mainfilename);
7511 if ( tPDB == NULL ) {
7512 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7516 NAMD_die(
"Number of atoms in rotating drag tag PDB doesn't match coordinate PDB");
7524 if (rotDragAxisFile == NULL) {
7525 if (rotDragFile == NULL) {
7526 NAMD_die(
"Rotating drag axis file can not be same as coordinate PDB file");
7528 if (rotDragAxisFile->
next != NULL) {
7529 NAMD_die(
"Multiple definitions of rotating drag axis file in configuration file");
7531 if (rotDragPivotFile == NULL) {
7532 NAMD_die(
"Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
7539 if ( (cwd == NULL) || (rotDragAxisFile->
data[0] ==
'/') ) {
7540 strcpy(axisfilename, rotDragAxisFile->
data);
7542 strcpy(axisfilename, cwd);
7543 strcat(axisfilename, rotDragAxisFile->
data);
7546 aPDB =
new PDB(axisfilename);
7547 if ( aPDB == NULL ) {
7548 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7552 NAMD_die(
"Number of atoms in rotating drag axis PDB doesn't match coordinate PDB");
7560 if (rotDragPivotFile == NULL) {
7561 if (rotDragFile == NULL) {
7562 NAMD_die(
"Rotating drag pivot point file can not be same as coordinate PDB file");
7564 if (rotDragPivotFile->
next != NULL) {
7565 NAMD_die(
"Multiple definitions of rotating drag pivot point file in configuration file");
7567 if (rotDragAxisFile == NULL) {
7568 NAMD_die(
"Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
7575 if ( (cwd == NULL) || (rotDragPivotFile->
data[0] ==
'/') ) {
7576 strcpy(pivotfilename, rotDragPivotFile->
data);
7578 strcpy(pivotfilename, cwd);
7579 strcat(pivotfilename, rotDragPivotFile->
data);
7582 pPDB =
new PDB(pivotfilename);
7583 if ( pPDB == NULL ) {
7584 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7588 NAMD_die(
"Number of atoms in rotating drag pivot point PDB doesn't match coordinate PDB");
7597 if (rotDragVelFile == NULL) {
7600 if (rotDragVelFile->
next != NULL) {
7601 NAMD_die(
"Multiple definitions of rotating drag velocity file in configuration file");
7604 if ( (cwd == NULL) || (rotDragVelFile->
data[0] ==
'/') ) {
7605 strcpy(velfilename, rotDragVelFile->
data);
7607 strcpy(velfilename, cwd);
7608 strcat(velfilename, rotDragVelFile->
data);
7611 vPDB =
new PDB(velfilename);
7612 if ( vPDB == NULL ) {
7613 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7617 NAMD_die(
"Number of atoms in rotating drag velocity PDB doesn't match coordinate PDB");
7628 if (rotDragCol == NULL) {
7631 if (rotDragCol->
next != NULL) {
7632 NAMD_die(
"Multiple definitions of drag tag column in config file");
7635 if ( rotDragFile == NULL
7636 && (!strcasecmp(rotDragCol->
data,
"X")
7637 || !strcasecmp(rotDragCol->
data,
"Y")
7638 || !strcasecmp(rotDragCol->
data,
"Z"))) {
7639 NAMD_die(
"Can not read rotating drag tags from X, Y, or Z column of the PDB coordinate file");
7641 if (!strcasecmp(rotDragCol->
data,
"X")) {
7643 }
else if (!strcasecmp(rotDragCol->
data,
"Y")) {
7645 }
else if (!strcasecmp(rotDragCol->
data,
"Z")) {
7647 }
else if (!strcasecmp(rotDragCol->
data,
"O")) {
7649 }
else if (!strcasecmp(rotDragCol->
data,
"B")) {
7653 NAMD_die(
"rotDragCol must have value of X, Y, Z, O, or B");
7665 if (rotDragVelCol == NULL) {
7668 if (rotDragVelCol->
next != NULL) {
7669 NAMD_die(
"Multiple definitions of drag angular velocity column in config file");
7672 if (rotDragVelFile == NULL
7673 && rotDragFile == NULL
7674 && strcasecmp(rotDragCol->
data,
"B")
7675 && strcasecmp(rotDragCol->
data,
"O")) {
7676 NAMD_die(
"Can not read rotating drag angular velocities from X, Y, or Z column of the PDB coordinate file");
7678 if (!strcasecmp(rotDragVelCol->
data,
"X")) {
7680 }
else if (!strcasecmp(rotDragVelCol->
data,
"Y")) {
7682 }
else if (!strcasecmp(rotDragVelCol->
data,
"Z")) {
7684 }
else if (!strcasecmp(rotDragVelCol->
data,
"O")) {
7686 }
else if (!strcasecmp(rotDragVelCol->
data,
"B")) {
7690 NAMD_die(
"rotDragVelCol must have value of X, Y, Z, O, or B");
7697 rotDragIndexes =
new int32[numAtoms];
7698 if (rotDragIndexes == NULL) {
7699 NAMD_die(
"memory allocation failed in Molecule::build_rotdrag_params()");
7703 for (i=0; i<numAtoms; i++) {
7706 dtval = (tPDB->
atom(i))->xcoor();
7709 dtval = (tPDB->
atom(i))->ycoor();
7712 dtval = (tPDB->
atom(i))->zcoor();
7715 dtval = (tPDB->
atom(i))->occupancy();
7718 dtval = (tPDB->
atom(i))->temperaturefactor();
7724 rotDragIndexes[i] = current_index;
7728 rotDragIndexes[i] = -1;
7732 if (current_index == 0) {
7733 iout <<
iWARN <<
"NO DRAGGED ATOMS WERE FOUND, BUT ROTATING DRAG IS ON . . . " <<
endi;
7735 rotDragParams =
new RotDragParams[current_index];
7736 if (rotDragParams == NULL) {
7737 NAMD_die(
"memory allocation failed in Molecule::build_rotdrag_params");
7741 numRotDrag = current_index;
7745 for (i=0; i<numAtoms; i++) {
7746 if (rotDragIndexes[i] != -1) {
7747 rotDragParams[rotDragIndexes[i]].a[0] = (aPDB->
atom(i))->xcoor();
7748 rotDragParams[rotDragIndexes[i]].a[1] = (aPDB->
atom(i))->ycoor();
7749 rotDragParams[rotDragIndexes[i]].a[2] = (aPDB->
atom(i))->zcoor();
7750 rotDragParams[rotDragIndexes[i]].p[0] = (pPDB->
atom(i))->xcoor();
7751 rotDragParams[rotDragIndexes[i]].p[1] = (pPDB->
atom(i))->ycoor();
7752 rotDragParams[rotDragIndexes[i]].p[2] = (pPDB->
atom(i))->zcoor();
7755 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->xcoor();
7758 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->ycoor();
7761 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->zcoor();
7764 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->occupancy();
7767 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->temperaturefactor();
7773 if (rotDragFile != NULL)
delete tPDB;
7774 if (rotDragAxisFile != NULL)
delete aPDB;
7775 if (rotDragPivotFile != NULL)
delete pPDB;
7776 if (rotDragVelFile != NULL)
delete vPDB;
7813 PDB *tPDB, *aPDB, *pPDB, *vPDB;
7815 int current_index=0;
7828 if (consTorqueFile == NULL) {
7829 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, consTorqueFile required.");
7834 if (consTorqueFile->
next != NULL) {
7835 NAMD_die(
"Multiple definitions of \"constant\" torque tag file in configuration file");
7838 if ( (cwd == NULL) || (consTorqueFile->
data[0] ==
'/') ) {
7839 strcpy(mainfilename, consTorqueFile->
data);
7841 strcpy(mainfilename, cwd);
7842 strcat(mainfilename, consTorqueFile->
data);
7845 tPDB =
new PDB(mainfilename);
7846 if ( tPDB == NULL ) {
7847 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7851 NAMD_die(
"Number of atoms in \"constant\" torque tag PDB doesn't match coordinate PDB");
7859 if (consTorqueAxisFile == NULL) {
7860 if (consTorqueFile == NULL) {
7861 NAMD_die(
"\"Constant\" torque axis file can not be same as coordinate PDB file");
7863 if (consTorqueAxisFile->
next != NULL) {
7864 NAMD_die(
"Multiple definitions of \"constant\" torque axis file in configuration file");
7866 if (consTorquePivotFile == NULL) {
7867 NAMD_die(
"Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
7874 if ( (cwd == NULL) || (consTorqueAxisFile->
data[0] ==
'/') ) {
7875 strcpy(axisfilename, consTorqueAxisFile->
data);
7877 strcpy(axisfilename, cwd);
7878 strcat(axisfilename, consTorqueAxisFile->
data);
7881 aPDB =
new PDB(axisfilename);
7882 if ( aPDB == NULL ) {
7883 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7887 NAMD_die(
"Number of atoms in \"constant\" torque axis PDB doesn't match coordinate PDB");
7895 if (consTorquePivotFile == NULL) {
7896 if (consTorqueFile == NULL) {
7897 NAMD_die(
"\"Constant\" torque pivot point file can not be same as coordinate PDB file");
7899 if (consTorquePivotFile->
next != NULL) {
7900 NAMD_die(
"Multiple definitions of \"constant\" torque pivot point file in configuration file");
7902 if (consTorqueAxisFile == NULL) {
7903 NAMD_die(
"Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
7910 if ( (cwd == NULL) || (consTorquePivotFile->
data[0] ==
'/') ) {
7911 strcpy(pivotfilename, consTorquePivotFile->
data);
7913 strcpy(pivotfilename, cwd);
7914 strcat(pivotfilename, consTorquePivotFile->
data);
7917 pPDB =
new PDB(pivotfilename);
7918 if ( pPDB == NULL ) {
7919 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7923 NAMD_die(
"Number of atoms in \"constant\" torque pivot point PDB doesn't match coordinate PDB");
7932 if (consTorqueValFile == NULL) {
7935 if (consTorqueValFile->
next != NULL) {
7936 NAMD_die(
"Multiple definitions of \"constant\" torque velocity file in configuration file");
7939 if ( (cwd == NULL) || (consTorqueValFile->
data[0] ==
'/') ) {
7940 strcpy(velfilename, consTorqueValFile->
data);
7942 strcpy(velfilename, cwd);
7943 strcat(velfilename, consTorqueValFile->
data);
7946 vPDB =
new PDB(velfilename);
7947 if ( vPDB == NULL ) {
7948 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7952 NAMD_die(
"Number of atoms in \"constant\" torque velocity PDB doesn't match coordinate PDB");
7963 if (consTorqueCol == NULL) {
7966 if (consTorqueCol->
next != NULL) {
7967 NAMD_die(
"Multiple definitions of torque tag column in config file");
7970 if ( consTorqueFile == NULL
7971 && (!strcasecmp(consTorqueCol->
data,
"X")
7972 || !strcasecmp(consTorqueCol->
data,
"Y")
7973 || !strcasecmp(consTorqueCol->
data,
"Z"))) {
7974 NAMD_die(
"Can not read \"constant\" torque tags from X, Y, or Z column of the PDB coordinate file");
7976 if (!strcasecmp(consTorqueCol->
data,
"X")) {
7978 }
else if (!strcasecmp(consTorqueCol->
data,
"Y")) {
7980 }
else if (!strcasecmp(consTorqueCol->
data,
"Z")) {
7982 }
else if (!strcasecmp(consTorqueCol->
data,
"O")) {
7984 }
else if (!strcasecmp(consTorqueCol->
data,
"B")) {
7988 NAMD_die(
"consTorqueCol must have value of X, Y, Z, O, or B");
8000 if (consTorqueValCol == NULL) {
8003 if (consTorqueValCol->
next != NULL) {
8004 NAMD_die(
"Multiple definitions of torque value column in config file");
8007 if (consTorqueValFile == NULL
8008 && consTorqueFile == NULL
8009 && strcasecmp(consTorqueCol->
data,
"B")
8010 && strcasecmp(consTorqueCol->
data,
"O")) {
8011 NAMD_die(
"Can not read \"constant\" torque values from X, Y, or Z column of the PDB coordinate file");
8013 if (!strcasecmp(consTorqueValCol->
data,
"X")) {
8015 }
else if (!strcasecmp(consTorqueValCol->
data,
"Y")) {
8017 }
else if (!strcasecmp(consTorqueValCol->
data,
"Z")) {
8019 }
else if (!strcasecmp(consTorqueValCol->
data,
"O")) {
8021 }
else if (!strcasecmp(consTorqueValCol->
data,
"B")) {
8025 NAMD_die(
"consTorqueValCol must have value of X, Y, Z, O, or B");
8032 consTorqueIndexes =
new int32[numAtoms];
8033 if (consTorqueIndexes == NULL) {
8034 NAMD_die(
"memory allocation failed in Molecule::build_constorque_params()");
8038 for (i=0; i<numAtoms; i++) {
8041 dtval = (tPDB->
atom(i))->xcoor();
8044 dtval = (tPDB->
atom(i))->ycoor();
8047 dtval = (tPDB->
atom(i))->zcoor();
8050 dtval = (tPDB->
atom(i))->occupancy();
8053 dtval = (tPDB->
atom(i))->temperaturefactor();
8059 consTorqueIndexes[i] = current_index;
8063 consTorqueIndexes[i] = -1;
8067 if (current_index == 0) {
8068 iout <<
iWARN <<
"NO TORQUED ATOMS WERE FOUND, BUT \"CONSTANT\" TORQUE IS ON . . . " <<
endi;
8070 consTorqueParams =
new ConsTorqueParams[current_index];
8071 if (consTorqueParams == NULL) {
8072 NAMD_die(
"memory allocation failed in Molecule::build_constorque_params");
8076 numConsTorque = current_index;
8080 for (i=0; i<numAtoms; i++) {
8081 if (consTorqueIndexes[i] != -1) {
8082 consTorqueParams[consTorqueIndexes[i]].a[0] = (aPDB->
atom(i))->xcoor();
8083 consTorqueParams[consTorqueIndexes[i]].a[1] = (aPDB->
atom(i))->ycoor();
8084 consTorqueParams[consTorqueIndexes[i]].a[2] = (aPDB->
atom(i))->zcoor();
8085 consTorqueParams[consTorqueIndexes[i]].p[0] = (pPDB->
atom(i))->xcoor();
8086 consTorqueParams[consTorqueIndexes[i]].p[1] = (pPDB->
atom(i))->ycoor();
8087 consTorqueParams[consTorqueIndexes[i]].p[2] = (pPDB->
atom(i))->zcoor();
8090 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->xcoor();
8093 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->ycoor();
8096 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->zcoor();
8099 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->occupancy();
8102 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->temperaturefactor();
8108 if (consTorqueFile != NULL)
delete tPDB;
8109 if (consTorqueAxisFile != NULL)
delete aPDB;
8110 if (consTorquePivotFile != NULL)
delete pPDB;
8111 if (consTorqueValFile != NULL)
delete vPDB;
8137 iout <<
iWARN <<
"NO CONSTANT FORCES SPECIFIED, BUT CONSTANT FORCE IS ON . . .\n" <<
endi;
8138 consForceIndexes =
new int32[numAtoms];
8139 for (i=0; i<numAtoms; i++) consForceIndexes[i] = -1;
8143 if ((forcePDB=
new PDB(filename)) == NULL)
8144 NAMD_die(
"Memory allocation failed in Molecule::build_constant_forces");
8146 NAMD_die(
"Number of atoms in constant force PDB doesn't match coordinate PDB");
8151 consForceIndexes =
new int32[numAtoms];
8152 if (consForceIndexes == NULL)
8153 NAMD_die(
"memory allocation failed in Molecule::build_constant_forces()");
8157 for (i=0; i<numAtoms; i++)
8161 consForceIndexes[i] = -1;
8164 consForceIndexes[i] = numConsForce++;
8166 if (numConsForce == 0)
8168 iout <<
iWARN <<
"NO NON-ZERO FORCES WERE FOUND, BUT CONSTANT FORCE IS ON . . .\n" <<
endi;
8171 consForce =
new Vector[numConsForce];
8172 if (consForce == NULL)
8173 NAMD_die(
"memory allocation failed in Molecule::build_constant_forces");
8175 for (i=0; i<numAtoms; i++)
8176 if ((index=consForceIndexes[i]) != -1)
8193 langevinParams =
new Real[numAtoms];
8195 if ( (langevinParams == NULL) )
8197 NAMD_die(
"memory allocation failed in Molecule::build_langevin_params()");
8201 for (
int i=0; i<numAtoms; i++)
8205 if ( (! doHydrogen) && is_hydrogen(i) ) bval = 0;
8206 else if ( is_lp(i) ) bval = 0;
8207 else if ( is_drude(i) ) bval = drudeCoupling;
8210 langevinParams[i] = bval;
8247 if (langfile == NULL)
8249 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, langevinFile required.");
8254 if (langfile->
next != NULL)
8256 NAMD_die(
"Multiple definitions of langvein PDB file in configuration file");
8259 if ( (cwd == NULL) || (langfile->
data[0] ==
'/') )
8261 strcpy(filename, langfile->
data);
8265 strcpy(filename, cwd);
8266 strcat(filename, langfile->
data);
8269 bPDB =
new PDB(filename);
8272 NAMD_die(
"Memory allocation failed in Molecule::build_langevin_params");
8277 NAMD_die(
"Number of atoms in langevin parameter PDB doesn't match coordinate PDB");
8286 if (langcol == NULL)
8292 if (langcol->
next != NULL)
8294 NAMD_die(
"Multiple definitions of langevin parameter column in config file");
8297 if (strcasecmp(langcol->
data,
"X") == 0)
8301 else if (strcasecmp(langcol->
data,
"Y") == 0)
8305 else if (strcasecmp(langcol->
data,
"Z") == 0)
8309 else if (strcasecmp(langcol->
data,
"O") == 0)
8313 else if (strcasecmp(langcol->
data,
"B") == 0)
8319 NAMD_die(
"langevincol must have value of X, Y, Z, O, or B");
8324 langevinParams =
new Real[numAtoms];
8326 if ( (langevinParams == NULL) )
8328 NAMD_die(
"memory allocation failed in Molecule::build_langevin_params()");
8332 for (i=0; i<numAtoms; i++)
8338 bval = (bPDB->
atom(i))->xcoor();
8341 bval = (bPDB->
atom(i))->ycoor();
8344 bval = (bPDB->
atom(i))->zcoor();
8347 bval = (bPDB->
atom(i))->occupancy();
8350 bval = (bPDB->
atom(i))->temperaturefactor();
8355 langevinParams[i] = bval;
8359 if (langfile != NULL)
8398 if (fixedfile == NULL)
8400 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, fixedAtomsFile required.");
8405 if (fixedfile->
next != NULL)
8407 NAMD_die(
"Multiple definitions of fixed atoms PDB file in configuration file");
8410 if ( (cwd == NULL) || (fixedfile->
data[0] ==
'/') )
8412 strcpy(filename, fixedfile->
data);
8416 strcpy(filename, cwd);
8417 strcat(filename, fixedfile->
data);
8420 bPDB =
new PDB(filename);
8423 NAMD_die(
"Memory allocation failed in Molecule::build_fixed_atoms");
8428 NAMD_die(
"Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
8437 if (fixedcol == NULL)
8443 if (fixedcol->
next != NULL)
8445 NAMD_die(
"Multiple definitions of fixed atoms column in config file");
8448 if (strcasecmp(fixedcol->
data,
"X") == 0)
8452 else if (strcasecmp(fixedcol->
data,
"Y") == 0)
8456 else if (strcasecmp(fixedcol->
data,
"Z") == 0)
8460 else if (strcasecmp(fixedcol->
data,
"O") == 0)
8464 else if (strcasecmp(fixedcol->
data,
"B") == 0)
8470 NAMD_die(
"fixedatomscol must have value of X, Y, Z, O, or B");
8475 fixedAtomFlags =
new int32[numAtoms];
8477 if (fixedAtomFlags == NULL)
8479 NAMD_die(
"memory allocation failed in Molecule::build_fixed_atoms()");
8485 for (i=0; i<numAtoms; i++)
8491 bval = (bPDB->
atom(i))->xcoor();
8494 bval = (bPDB->
atom(i))->ycoor();
8497 bval = (bPDB->
atom(i))->zcoor();
8500 bval = (bPDB->
atom(i))->occupancy();
8503 bval = (bPDB->
atom(i))->temperaturefactor();
8509 fixedAtomFlags[i] = 1;
8513 fixedAtomFlags[i] = 0;
8518 if (fixedfile != NULL)
8525 if ( numRigidBonds ) {
8527 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
8528 int parentIsFixed = 0;
8529 for( ; h_i != h_e; ++h_i ) {
8531 parentIsFixed = fixedAtomFlags[h_i->
atomID];
8532 if ( (rigidBondLengths[h_i->
atomID] > 0.)
8533 && fixedAtomFlags[h_i[1].atomID]
8534 && fixedAtomFlags[h_i[2].atomID] ) {
8535 ++numFixedRigidBonds;
8538 if ( (rigidBondLengths[h_i->
atomID] > 0.)
8539 && fixedAtomFlags[h_i->
atomID]
8540 && parentIsFixed ) {
8541 ++numFixedRigidBonds;
8551 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
8553 for( ; h_i != h_e; ++h_i ) {
8555 if ( allFixed ) ++numFixedGroups;
8558 allFixed = allFixed && fixedAtomFlags[h_i->
atomID];
8560 if ( allFixed ) ++numFixedGroups;
8597 int current_index=0;
8605 if (stirredfile == NULL)
8607 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, stirFilename required.");
8610 iout <<
iWARN <<
"STIRRING USING INITIAL POSITION FILE FOR REFERENCE POSITIONS" <<
endi;
8614 if (stirredfile->
next != NULL)
8616 NAMD_die(
"Multiple definitions of stirred atoms PDB file in configuration file");
8619 if ( (cwd == NULL) || (stirredfile->
data[0] ==
'/') )
8621 strcpy(filename, stirredfile->
data);
8625 strcpy(filename, cwd);
8626 strcat(filename, stirredfile->
data);
8630 sPDB =
new PDB(filename);
8634 NAMD_die(
"Memory allocation failed in Molecule::build_stirred_atoms");
8640 NAMD_die(
"Number of atoms in stirred atoms PDB doesn't match coordinate PDB");
8652 if (stirredcol == NULL)
8658 if (stirredcol->
next != NULL)
8660 NAMD_die(
"Multiple definitions of stirred atoms column in config file");
8663 if (strcasecmp(stirredcol->
data,
"O") == 0)
8667 else if (strcasecmp(stirredcol->
data,
"B") == 0)
8673 NAMD_die(
"stirredAtomsCol must have value of O or B");
8680 stirIndexes =
new int32[numAtoms];
8682 if (stirIndexes == NULL)
8684 NAMD_die(
"memory allocation failed in Molecule::build_stirred_params()");
8689 for (i=0; i<numAtoms; i++)
8699 bval = (sPDB->
atom(i))->occupancy();
8702 bval = (sPDB->
atom(i))->temperaturefactor();
8711 stirIndexes[i] = current_index;
8717 stirIndexes[i] = -1;
8725 if (current_index == 0)
8728 iout <<
iWARN <<
"NO STIRRED ATOMS WERE FOUND, BUT STIRRING TORQUES ARE ON . . .\n" <<
endi;
8733 stirParams =
new StirParams[current_index];
8735 if (stirParams == NULL)
8737 NAMD_die(
"memory allocation failed in Molecule::build_stir_params");
8741 numStirredAtoms = current_index;
8745 for (i=0; i<numAtoms; i++)
8747 if (stirIndexes[i] != -1)
8751 stirParams[stirIndexes[i]].refPos.x = (sPDB->
atom(i))->xcoor();
8752 stirParams[stirIndexes[i]].refPos.y = (sPDB->
atom(i))->ycoor();
8753 stirParams[stirIndexes[i]].refPos.z = (sPDB->
atom(i))->zcoor();
8758 if (stirredfile != NULL)
8774 int a1,a2,a3,a4;
float k, ref, upper;
8775 int anglesNormal = (
simParams->extraBondsCosAngles ? 0 : 1 );
8776 #ifndef MEM_OPT_VERSION 8789 NAMD_die(
"NO EXTRA BONDS FILES SPECIFIED");
8792 for ( ; file; file = file->
next ) {
8793 FILE *f = fopen(file->
data,
"r");
8795 sprintf(err_msg,
"UNABLE TO OPEN EXTRA BONDS FILE %s", file->
data);
8807 if (ret_code!=0)
break;
8810 sscanf(buffer,
"%s",type);
8812 #define CHECKATOMID(ATOMID) if ( ATOMID < 0 || ATOMID >= numAtoms ) badatom = 1; 8816 if ( ! strncasecmp(type,
"bond",4) ) {
8817 if ( sscanf(buffer,
"%s %d %d %f %f %s",
8818 type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1;
8824 #ifndef MEM_OPT_VERSION 8832 tmpv.
k = k; tmpv.
x0 = ref;
8834 }
else if ( ! strncasecmp(type,
"wall",4) ) {
8837 if ( sscanf(buffer,
"%s %d %d %f %f %f %s",
8838 type, &a1, &a2, &k, &ref, &upper, err_msg) != 6 ) badline = 1;
8839 else if (upper < ref) badline = 1;
8845 #ifndef MEM_OPT_VERSION 8853 tmpv.
k = k; tmpv.
x0 = ref; tmpv.
x1 = upper;
8855 }
else if ( ! strncasecmp(type,
"angle",4) ) {
8856 if ( sscanf(buffer,
"%s %d %d %d %f %f %s",
8857 type, &a1, &a2, &a3, &k, &ref, err_msg) != 6 ) badline = 1;
8863 #ifndef MEM_OPT_VERSION 8871 tmpv.
k = k; tmpv.
theta0 = ref / 180. *
PI;
8873 tmpv.
normal = anglesNormal;
8876 }
else if ( ! strncasecmp(type,
"dihedral",4) ) {
8878 int ret = 1 + sscanf(buffer,
"%s %d %d %d %d %f %f %s",
8879 type, &a1, &a2, &a3, &a4, &k, &ref, err_msg);
8881 ret = sscanf(buffer,
"%s %d %d %d %d %f %d %f %s",
8882 type, &a1, &a2, &a3, &a4, &k, &n, &ref, err_msg);
8884 if ( ret != 8 ) badline = 1;
8891 #ifndef MEM_OPT_VERSION 8902 }
else if ( ! strncasecmp(type,
"improper",4) ) {
8904 int ret = 1 + sscanf(buffer,
"%s %d %d %d %d %f %f %s",
8905 type, &a1, &a2, &a3, &a4, &k, &ref, err_msg);
8907 ret = sscanf(buffer,
"%s %d %d %d %d %f %d %f %s",
8908 type, &a1, &a2, &a3, &a4, &k, &n, &ref, err_msg);
8910 if ( ret != 8 ) badline = 1;
8917 #ifndef MEM_OPT_VERSION 8928 }
else if ( ! strncasecmp(type,
"#",1) ) {
8935 sprintf(err_msg,
"BAD LINE IN EXTRA BONDS FILE %s: %s",
8936 file->
data, buffer);
8940 sprintf(err_msg,
"BAD ATOM ID IN EXTRA BONDS FILE %s: %s",
8941 file->
data, buffer);
8950 if ( extraNumBonds ) {
8951 iout <<
iINFO <<
"READ " << extraNumBonds <<
" EXTRA BONDS\n" <<
endi;
8953 #ifndef MEM_OPT_VERSION 8954 Bond *newbonds =
new Bond[numBonds+extraNumBonds];
8955 memcpy(newbonds, this->bonds, numBonds*
sizeof(
Bond));
8956 memcpy(newbonds+numBonds, bonds.
begin(), extraNumBonds*
sizeof(
Bond));
8957 delete [] this->bonds;
8958 this->bonds = newbonds;
8959 numBonds += extraNumBonds;
8974 if ( extraNumAngles ) {
8975 iout <<
iINFO <<
"READ " << extraNumAngles <<
" EXTRA ANGLES\n" <<
endi;
8976 if ( anglesNormal ) {
8977 iout <<
iINFO <<
"EXTRA ANGLES ARE NORMAL HARMONIC\n" <<
endi;
8978 }
else if (
simParams->extraBondsCosAnglesSetByUser ) {
8981 iout <<
iWARN <<
"EXTRA ANGLES ARE COSINE-BASED BY DEFAULT TO MATCH PREVIOUS VERSIONS\n";
8982 iout <<
iWARN <<
"FOR NORMAL HARMONIC EXTRA ANGLES SET extraBondsCosAngles off\n" <<
endi;
8984 #ifndef MEM_OPT_VERSION 8985 Angle *newangles =
new Angle[numAngles+extraNumAngles];
8986 memcpy(newangles, this->angles, numAngles*
sizeof(
Angle));
8987 memcpy(newangles+numAngles, angles.
begin(), extraNumAngles*
sizeof(
Angle));
8988 delete [] this->angles;
8989 this->angles = newangles;
8990 numAngles += extraNumAngles;
9005 if ( extraNumDihedrals ) {
9006 iout <<
iINFO <<
"READ " << extraNumDihedrals <<
" EXTRA DIHEDRALS\n" <<
endi;
9007 #ifndef MEM_OPT_VERSION 9009 memcpy(newdihedrals, this->dihedrals, numDihedrals*
sizeof(
Dihedral));
9010 memcpy(newdihedrals+numDihedrals, dihedrals.
begin(), extraNumDihedrals*
sizeof(
Dihedral));
9011 delete [] this->dihedrals;
9012 this->dihedrals = newdihedrals;
9013 numDihedrals += extraNumDihedrals;
9028 if ( extraNumImpropers ) {
9029 iout <<
iINFO <<
"READ " << extraNumImpropers <<
" EXTRA IMPROPERS\n" <<
endi;
9030 #ifndef MEM_OPT_VERSION 9032 memcpy(newimpropers, this->impropers, numImpropers*
sizeof(
Improper));
9033 memcpy(newimpropers+numImpropers, impropers.
begin(), extraNumImpropers*
sizeof(
Improper));
9034 delete [] this->impropers;
9035 this->impropers = newimpropers;
9036 numImpropers += extraNumImpropers;
9070 PDB *initial_pdb,
char *cwd,
9071 const char *simmethod) {
9081 if (alchfile == NULL) {
9082 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, alchfile required.");
9084 strcpy(filename,
"coordinate pdb file (default)");
9087 if (alchfile->
next != NULL) {
9088 char *new_err_msg =
new char[24 + strlen(simmethod) + 26];
9089 sprintf(new_err_msg,
"Multiple definitions of %sFile in configuration file",simmethod);
9093 if ((cwd == NULL) || (alchfile->
data[0] ==
'/')) {
9094 strcpy(filename, alchfile->
data);
9097 strcpy(filename, cwd);
9098 strcat(filename, alchfile->
data);
9101 bPDB =
new PDB(filename);
9103 NAMD_die(
"Memory allocation failed in Molecule:build_fep_flags");
9107 char *new_err_msg =
new char[19 + strlen(simmethod) + 38];
9108 sprintf(new_err_msg,
"Number of atoms in %sFile PDB does not match coordinate PDB",simmethod);
9116 if (alchcol == NULL) {
9120 if (alchcol->
next != NULL) {
9121 char *new_err_msg =
new char[24 + strlen(simmethod) + 35];
9122 sprintf(new_err_msg,
"Multiple definitions of %s parameter column in config file",simmethod);
9126 if (strcasecmp(alchcol->
data,
"X") == 0) {
9129 else if (strcasecmp(alchcol->
data,
"Y") == 0) {
9132 else if (strcasecmp(alchcol->
data,
"Z") == 0) {
9135 else if (strcasecmp(alchcol->
data,
"O") == 0) {
9138 else if (strcasecmp(alchcol->
data,
"B") == 0) {
9142 NAMD_die(
"alchcol must have value of X, Y, Z, O or B");
9146 iout <<
iINFO <<
"To read " << simmethod <<
"data from file: " << filename
9148 iout <<
iINFO <<
"To read " << simmethod <<
"flag data from column: " << bcol
9152 fepAtomFlags =
new unsigned char[numAtoms];
9154 if (fepAtomFlags == NULL) {
9155 NAMD_die(
"Memory allocation failed in Molecule::build_fep_params()");
9158 double lesMassFactor = 1.0;
9160 lesMassFactor = 1.0 /
simParams->lesFactor;
9164 for (i = 0; i < numAtoms; i++) {
9168 bval = (bPDB->
atom(i))->xcoor();
9171 bval = (bPDB->
atom(i))->ycoor();
9174 bval = (bPDB->
atom(i))->zcoor();
9177 bval = (bPDB->
atom(i))->occupancy();
9180 bval = (bPDB->
atom(i))->temperaturefactor();
9186 if ( bval == (
int) bval && bval > 0 ) {
9188 NAMD_die(
"LES flag must be less than or equal to lesFactor.");
9189 fepAtomFlags[i] = (int) bval;
9190 #ifdef MEM_OPT_VERSION 9191 Real newMass = atomMassPool[eachAtomMass[i]]*lesMassFactor;
9192 eachAtomMass[i] = insert_new_mass(newMass);
9194 atoms[i].mass *= lesMassFactor;
9199 fepAtomFlags[i] = 0;
9203 fepAtomFlags[i] = 3;
9205 }
else if (bval == 1.0) {
9208 }
else if (bval == -1.0) {
9209 fepAtomFlags[i] = 2;
9211 }
else if (bval == -2.0) {
9212 fepAtomFlags[i] = 4;
9215 fepAtomFlags[i] = 0;
9217 }
else if (
simParams->pairInteractionOn) {
9218 if (bval ==
simParams->pairInteractionGroup1) {
9219 fepAtomFlags[i] = 1;
9221 }
else if (bval ==
simParams->pairInteractionGroup2) {
9222 fepAtomFlags[i] = 2;
9225 fepAtomFlags[i] = 0;
9227 }
else if (
simParams->pressureProfileAtomTypes > 1) {
9228 fepAtomFlags[i] = (int) bval;
9230 #ifdef OPENATOM_VERSION 9234 fepAtomFlags[i] = bval;
9237 fepAtomFlags[i] = 0;
9240 #endif //OPENATOM_VERSION 9244 if (alchfile != NULL) {
9268 if (ssfile == NULL) {
9269 if ( ! initial_pdb ) {
9270 NAMD_die(
"Initial PDB file unavailable, soluteScalingFile required.");
9273 strcpy(filename,
"coordinate PDB file (default)");
9276 if (ssfile->
next != NULL) {
9277 NAMD_die(
"Multiple definitions of soluteScalingFile in configuration file");
9280 if ((cwd == NULL) || (ssfile->
data[0] ==
'/')) {
9281 strcpy(filename, ssfile->
data);
9284 strcpy(filename, cwd);
9285 strcat(filename, ssfile->
data);
9288 bPDB =
new PDB(filename);
9290 NAMD_die(
"Memory allocation failed in Molecule::build_ss_flags");
9294 NAMD_die(
"Number of atoms in soluteScalingFile PDB does not match coordinate PDB");
9298 if (sscol == NULL) {
9302 if (sscol->
next != NULL) {
9303 NAMD_die(
"Multiple definitions of soluteScalingCol value in config file");
9306 if (strcasecmp(sscol->
data,
"X") == 0) {
9309 else if (strcasecmp(sscol->
data,
"Y") == 0) {
9312 else if (strcasecmp(sscol->
data,
"Z") == 0) {
9315 else if (strcasecmp(sscol->
data,
"O") == 0) {
9318 else if (strcasecmp(sscol->
data,
"B") == 0) {
9322 NAMD_die(
"soluteScalingCol must have value of X, Y, Z, O or B");
9326 iout <<
iINFO <<
"Reading solute scaling data from file: " 9327 << filename <<
"\n" <<
endi;
9328 iout <<
iINFO <<
"Reading solute scaling flags from column: " 9329 << bcol <<
"\n" <<
endi;
9331 ssAtomFlags =
new unsigned char[numAtoms];
9332 ss_index =
new int[numAtoms];
9334 if (ssAtomFlags == NULL || ss_index == NULL) {
9335 NAMD_die(
"Memory allocation failed in Molecule::build_ss_params()");
9339 for (i = 0; i < numAtoms; i++) {
9342 bval = (bPDB->
atom(i))->xcoor();
9345 bval = (bPDB->
atom(i))->ycoor();
9348 bval = (bPDB->
atom(i))->zcoor();
9351 bval = (bPDB->
atom(i))->occupancy();
9354 bval = (bPDB->
atom(i))->temperaturefactor();
9360 ss_index[num_ss] = i;
9369 if (ssfile != NULL) {
9378 int *numAtomsByLjType =
new int[LJtypecount];
9382 ss_vdw_type =
new int[LJtypecount];
9385 for (i = 0; i < LJtypecount; i++) {
9386 numAtomsByLjType[i] = 0;
9391 for (i = 0; i < num_ss; i++) {
9392 numAtomsByLjType[atoms[ss_index[i]].vdw_type]++;
9396 ss_num_vdw_params = 0;
9397 for (i = 0; i < LJtypecount; i++) {
9399 if (numAtomsByLjType[i] != 0) {
9402 ss_vdw_type[ss_num_vdw_params] = i;
9405 ss_num_vdw_params++;
9409 for (i = 0; i < num_ss; i++) {
9411 for (j = 0; j < ss_num_vdw_params; j++) {
9413 if (atoms[ss_index[i]].vdw_type == ss_vdw_type[j]) {
9416 atoms[ss_index[i]].vdw_type = LJtypecount + j;
9421 delete [] numAtomsByLjType;
9435 #ifndef MEM_OPT_VERSION 9439 FILE *alch_unpert_bond_file;
9442 if ((alch_unpert_bond_file =
Fopen(alch_fname,
"r")) == NULL) {
9443 sprintf(err_msg,
"UNABLE TO OPEN ALCH UNPERTBURBED BOND FILE %s", alch_fname);
9453 NAMD_die(
"DID NOT FIND NBOND AFTER ATOM LIST IN ALCH UNPERT PSF");
9457 sscanf(buffer,
"%d", &num_alch_unpert_Bonds);
9459 read_alch_unpert_bonds(alch_unpert_bond_file);
9468 NAMD_die(
"DID NOT FIND NTHETA AFTER BOND LIST IN ALCH UNPERT PSF");
9472 sscanf(buffer,
"%d", &num_alch_unpert_Angles);
9474 read_alch_unpert_angles(alch_unpert_bond_file);
9485 NAMD_die(
"DID NOT FIND NPHI AFTER ANGLE LIST IN ALCH UNPERT PSF");
9489 sscanf(buffer,
"%d", &num_alch_unpert_Dihedrals);
9491 read_alch_unpert_dihedrals(alch_unpert_bond_file);
9493 Fclose(alch_unpert_bond_file);
9499 suspiciousAlchBonds = 0;
9500 for (
int i = 0; i < numBonds; i++) {
9501 int part1 = fepAtomFlags[bonds[i].atom1];
9502 int part2 = fepAtomFlags[bonds[i].atom2];
9503 if ((part1 == 1 || part2 == 1 ) &&
9504 (part1 == 2 || part2 == 2 )) {
9506 suspiciousAlchBonds++;
9511 Angle *nonalchAngles;
9512 nonalchAngles =
new Angle[numAngles];
9513 int nonalchAngleCount = 0;
9514 alchDroppedAngles = 0;
9515 for (
int i = 0; i < numAngles; i++) {
9516 int part1 = fepAtomFlags[angles[i].atom1];
9517 int part2 = fepAtomFlags[angles[i].atom2];
9518 int part3 = fepAtomFlags[angles[i].atom3];
9519 if ((part1 == 1 || part2 == 1 || part3 == 1) &&
9520 (part1 == 2 || part2 == 2 || part3 == 2)) {
9522 alchDroppedAngles++;
9525 if ( angles[i].angle_type == -1 ) {
9528 "MISSING PARAMETERS FOR ANGLE %i %i %i PARTITIONS %i %i %i\n",
9529 angles[i].atom1+1, angles[i].atom2+1, angles[i].atom3+1,
9530 part1, part2, part3);
9533 nonalchAngles[nonalchAngleCount++] = angles[i];
9536 numAngles = nonalchAngleCount;
9538 angles =
new Angle[numAngles];
9539 for (
int i = 0; i < nonalchAngleCount; i++) {
9540 angles[i]=nonalchAngles[i];
9542 delete [] nonalchAngles;
9547 nonalchDihedrals =
new Dihedral[numDihedrals];
9548 int nonalchDihedralCount = 0;
9549 alchDroppedDihedrals = 0;
9550 for (
int i = 0; i < numDihedrals; i++) {
9551 int part1 = fepAtomFlags[dihedrals[i].atom1];
9552 int part2 = fepAtomFlags[dihedrals[i].atom2];
9553 int part3 = fepAtomFlags[dihedrals[i].atom3];
9554 int part4 = fepAtomFlags[dihedrals[i].atom4];
9555 if ((part1 == 1 || part2 == 1 || part3 == 1 || part4 == 1) &&
9556 (part1 == 2 || part2 == 2 || part3 == 2 || part4 == 2)) {
9558 alchDroppedDihedrals++;
9561 if ( dihedrals[i].dihedral_type == -1 ) {
9564 "MISSING PARAMETERS FOR DIHEDRAL %i %i %i %i PARTITIONS %i %i %i %i\n",
9565 dihedrals[i].atom1+1, dihedrals[i].atom2+1,
9566 dihedrals[i].atom3+1, dihedrals[i].atom4+1,
9567 part1, part2, part3, part4);
9570 nonalchDihedrals[nonalchDihedralCount++] = dihedrals[i];
9573 numDihedrals = nonalchDihedralCount;
9574 delete [] dihedrals;
9575 dihedrals =
new Dihedral[numDihedrals];
9576 for (
int i = 0; i < numDihedrals; i++) {
9577 dihedrals[i]=nonalchDihedrals[i];
9579 delete [] nonalchDihedrals;
9583 nonalchImpropers =
new Improper[numImpropers];
9584 int nonalchImproperCount = 0;
9585 alchDroppedImpropers = 0;
9586 for (
int i = 0; i < numImpropers; i++) {
9587 int part1 = fepAtomFlags[impropers[i].atom1];
9588 int part2 = fepAtomFlags[impropers[i].atom2];
9589 int part3 = fepAtomFlags[impropers[i].atom3];
9590 int part4 = fepAtomFlags[impropers[i].atom4];
9591 if ((part1 == 1 || part2 == 1 || part3 == 1 || part4 == 1) &&
9592 (part1 == 2 || part2 == 2 || part3 == 2 || part4 == 2)) {
9594 alchDroppedImpropers++;
9597 nonalchImpropers[nonalchImproperCount++] = impropers[i];
9600 numImpropers = nonalchImproperCount;
9601 delete [] impropers;
9602 impropers =
new Improper[numImpropers];
9603 for (
int i = 0; i < numImpropers; i++) {
9604 impropers[i]=nonalchImpropers[i];
9606 delete [] nonalchImpropers;
9627 if (fixedfile == NULL) {
9628 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, excludeFromPressureFile required.");
9631 if (fixedfile->
next != NULL) {
9632 NAMD_die(
"Multiple definitions of excluded pressure atoms PDB file in configuration file");
9635 if ( (cwd == NULL) || (fixedfile->
data[0] ==
'/') ) {
9636 strcpy(filename, fixedfile->
data);
9638 strcpy(filename, cwd);
9639 strcat(filename, fixedfile->
data);
9641 bPDB =
new PDB(filename);
9642 if ( bPDB == NULL ) {
9643 NAMD_die(
"Memory allocation failed in Molecule::build_exPressure_atoms");
9647 NAMD_die(
"Number of atoms in excludedPressure atoms PDB doesn't match coordinate PDB");
9656 if (fixedcol == NULL) {
9659 if (fixedcol->
next != NULL) {
9660 NAMD_die(
"Multiple definitions of excludedPressure atoms column in config file");
9663 if (strcasecmp(fixedcol->
data,
"X") == 0) {
9665 }
else if (strcasecmp(fixedcol->
data,
"Y") == 0) {
9667 }
else if (strcasecmp(fixedcol->
data,
"Z") == 0) {
9669 }
else if (strcasecmp(fixedcol->
data,
"O") == 0) {
9671 }
else if (strcasecmp(fixedcol->
data,
"B") == 0) {
9674 NAMD_die(
"excludedPressureFileCol must have value of X, Y, Z, O, or B");
9679 exPressureAtomFlags =
new int32[numAtoms];
9681 if (exPressureAtomFlags == NULL) {
9682 NAMD_die(
"memory allocation failed in Molecule::build_fixed_atoms()");
9685 numExPressureAtoms = 0;
9688 for (i=0; i<numAtoms; i++) {
9691 case 1: bval = (bPDB->
atom(i))->xcoor();
break;
9692 case 2: bval = (bPDB->
atom(i))->ycoor();
break;
9693 case 3: bval = (bPDB->
atom(i))->zcoor();
break;
9694 case 4: bval = (bPDB->
atom(i))->occupancy();
break;
9695 case 5: bval = (bPDB->
atom(i))->temperaturefactor();
break;
9700 exPressureAtomFlags[i] = 1;
9701 numExPressureAtoms++;
9703 exPressureAtomFlags[i] = 0;
9706 if (fixedfile != NULL)
9709 iout <<
iINFO <<
"Got " << numExPressureAtoms <<
" excluded pressure atoms." 9719 return ((atoms[anum].status &
DrudeAtom) != 0);
9729 return ((atoms[anum].status &
OxygenAtom) != 0);
9734 return (hydrogenGroup[atoms[anum].hydrogenList].isGP);
9739 return (hydrogenGroup[atoms[anum].hydrogenList].waterVal == 2);
9744 return (hydrogenGroup[atoms[anum].hydrogenList].atomsInGroup);
9751 return atoms[anum].partner;
9755 if ( n != numAtoms )
9756 NAMD_die(
"Incorrect number of atoms in Molecule::reloadCharges().");
9758 #ifdef MEM_OPT_VERSION 9759 delete [] atomChargePool;
9760 vector<Real> tmpCharges;
9761 for(
int i=0; i<numAtoms; i++){
9765 for(
int j=0; j<tmpCharges.size();j++){
9766 if(tmpCharges[j] == charge[i]){
9772 tmpCharges.push_back(charge[i]);
9773 foundIdx = tmpCharges.size()-1;
9775 eachAtomCharge[i] = (
Index)foundIdx;
9777 chargePoolSize = tmpCharges.size();
9778 atomChargePool =
new Real[chargePoolSize];
9779 for(
int i=0; i<chargePoolSize; i++)
9780 atomChargePool[i] = tmpCharges[i];
9782 for(
int i=0; i<n; ++i ) atoms[i].charge = charge[i];
9794 #ifndef MEM_OPT_VERSION 9795 if(std::filesystem::path(dcdSelectionInputFile).extension() ==
".idx")
9798 IndexFile dcdSelectionInputIdx(dcdSelectionInputFile);
9799 std::vector<uint32> indexVec = dcdSelectionInputIdx.getAllElements();
9800 dcdSelectionParams[index].size = dcdSelectionInputIdx.size();
9801 int bitmask = 1 << index;
9802 for(
int i=0; i<indexVec.size(); ++i)
9803 { atoms[indexVec[i]].flags.dcdSelection |= bitmask; }
9807 PDB dcdSelectionInputFilePdb(dcdSelectionInputFile);
9809 int numDcdSelectionAtoms = dcdSelectionInputFilePdb.num_atoms();
9810 dcdSelectionParams[index].size = numDcdSelectionAtoms;
9811 if (numDcdSelectionAtoms < 1)
9812 NAMD_die(
"No atoms found in dcdSelectionInputFile\n");
9817 if (numDcdSelectionAtoms != numAtoms)
9818 NAMD_die(
"The number of atoms in dcdSelectionInputFile must be equal to the total number of atoms in the structure!");
9819 int bitmask = 1 << index;
9820 for (
int i=0; i<numAtoms; i++) {
9821 PDBAtom *atom = dcdSelectionInputFilePdb.atom(i);
9824 atoms[i].flags.dcdSelection |= bitmask;
9829 if(std::filesystem::path(dcdSelectionInputFile).extension() ==
".idx")
9832 IndexFile dcdSelectionInputIdx(dcdSelectionInputFile);
9833 dcdSelectionParams[index].size=dcdSelectionInputIdx.size();
9837 dcdSelectionParams[index].tag = index+1;
9844 #ifndef MEM_OPT_VERSION 9845 dcdSelectionParams[index].dcdSelectionIndexReverse.resize(numAtoms);
9847 int bitmask = 1 << index;
9849 for (
int i=0; i<numAtoms; i++) {
9850 if(atoms[i].flags.dcdSelection & bitmask)
9852 dcdSelectionParams[index].dcdSelectionIndex.push_back(i);
9853 dcdSelectionParams[index].dcdSelectionIndexReverse[i] = count++;
9857 dcdSelectionParams[index].dcdSelectionIndexReverse[i] = -1;
9860 dcdSelectionParams[index].size = count;
9868 dcdSelectionParams[index].frequency = freq;
9873 if(
auto keyIt=dcdSelectionKeyMap.find(std::string(keystr)); keyIt!=dcdSelectionKeyMap.end())
9875 return(keyIt->second);
9883 uint16 key=find_dcd_selection_index(keystr);
9890 key = dcdSelectionKeyMap.size();
9894 sprintf(errmsg,
"Key %s beyond Max (16) DCD Selections\n", keystr);
9897 dcdSelectionKeyMap[std::string(keystr)] = key;
9913 char *keystr =
new char[256];
9914 char *valstr =
new char[256];
9917 int dcdSelectionCounter[3]={0,0,0};
9918 StringList *dcdSelectionInputFileList = config->
find(
"dcdselectioninputfile");
9919 for(;dcdSelectionInputFileList; dcdSelectionInputFileList = dcdSelectionInputFileList->
next)
9921 sscanf(dcdSelectionInputFileList->
data,
"%80s %255s",keystr,valstr);
9922 key = find_or_create_dcd_selection_index(keystr);
9923 dcdSelectionCounter[0]++;
9924 build_dcd_selection_list_pdb(key, valstr);
9925 iout <<
iINFO <<
"Dcdselection tag "<< keystr <<
" Input file " << valstr <<
" index "<< key <<
"\n";
9927 StringList *dcdSelectionFileList = config->
find(
"dcdselectionfile");
9928 for(;dcdSelectionFileList; dcdSelectionFileList = dcdSelectionFileList->
next)
9931 sscanf(dcdSelectionFileList->
data,
"%80s %255s",keystr,valstr);
9932 key = find_dcd_selection_index(keystr);
9935 add_dcd_selection_file(key, valstr);
9936 dcdSelectionCounter[1]++;
9940 sprintf(errmsg,
"Keystring %s has no dcd input file match\n", keystr);
9943 iout <<
iINFO <<
"Dcdselection tag "<< keystr <<
" Output file " << valstr<<
" index "<< key <<
"\n";
9945 StringList *dcdSelectionFreqList = config->
find(
"dcdselectionfreq");
9946 for(;dcdSelectionFreqList; dcdSelectionFreqList = dcdSelectionFreqList->
next)
9949 sscanf(dcdSelectionFreqList->
data,
"%s %d",keystr,&freq);
9950 key = find_dcd_selection_index(keystr);
9953 add_dcd_selection_freq(key, freq);
9954 dcdSelectionCounter[2]++;
9958 sprintf(errmsg,
"Keystring %s has no dcd input file match\n", keystr);
9961 iout <<
iINFO <<
"Dcd Selection tag "<< keystr <<
" freq " << freq<<
" index "<< key <<
"\n";
9965 if((dcdSelectionCounter[0]!=dcdSelectionCounter[1] && CmiNumPartitions()==1) || dcdSelectionCounter[0]!=dcdSelectionCounter[2] || dcdSelectionCounter[0]>15)
9968 sprintf(errmsg,
"Invalid DCD Selection configuration\n");
9974 #ifndef MEM_OPT_VERSION 9981 void Molecule::build_atom_status(
void) {
9984 int numDrudeWaters = 0;
9987 numLonepairs = numZeroMassAtoms;
9989 iout <<
iWARN <<
"CORRECTION OF ZERO MASS ATOMS TURNED OFF " 9990 "BECAUSE LONE PAIRS ARE USED\n" <<
endi;
9994 if (is_lonepairs_psf && numLonepairs != numLphosts) {
9995 NAMD_die(
"must have same number of LP hosts as lone pairs");
9997 }
else if (numZeroMassAtoms) {
9998 for (i=0; i < numAtoms; i++) {
9999 if ( atoms[i].mass < 0.001 ) atoms[i].mass = 0.001;
10001 if ( ! CkMyPe() ) {
10002 iout <<
iWARN <<
"FOUND " << numZeroMassAtoms <<
10003 " ATOMS WITH ZERO OR NEGATIVE MASSES! CHANGED TO 0.001\n" <<
endi;
10008 hydrogenGroup.resize(numAtoms);
10010 for (i=0; i < numAtoms; i++) {
10011 atoms[i].partner = (-1);
10021 int hhbondcount = 0;
10022 for (i=0; i < numRealBonds; i++) {
10023 a1 = bonds[i].atom1;
10024 a2 = bonds[i].atom2;
10025 if (is_hydrogen(a1) && is_hydrogen(a2)) {
10028 atoms[a1].partner = a2;
10029 atoms[a2].partner = a1;
10037 if ( hhbondcount && ! CkMyPe() ) {
10038 iout <<
iWARN <<
"Found " << hhbondcount <<
" H-H bonds.\n" <<
endi;
10043 for (i=0; i < numRealBonds; i++) {
10044 a1 = bonds[i].atom1;
10045 a2 = bonds[i].atom2;
10046 if (is_hydrogen(a1)) {
10047 if (is_hydrogen(a2))
continue;
10048 atoms[a1].partner = a2;
10054 if (is_oxygen(a2)) hg[a2].waterVal++;
10056 if (is_hydrogen(a2)) {
10057 atoms[a2].partner = a1;
10063 if (is_oxygen(a1)) hg[a1].waterVal++;
10069 atoms[a1].partner = a2;
10076 atoms[a2].partner = a1;
10084 else if ( is_lonepairs_psf || is_drude_psf ) {
10085 if (is_lp(a1) || is_drude(a1)) {
10086 if (is_hydrogen(a2) || is_lp(a2) || is_drude(a2)) {
10088 sprintf(msg,
"%s particle %d is bonded to non-parent atom %d",
10089 (is_lp(a1) ?
"Lone pair" :
"Drude"), a1+1, a2+1);
10092 atoms[a1].partner = a2;
10098 else if (is_lp(a2) || is_drude(a2)) {
10099 if (is_hydrogen(a1) || is_lp(a1) || is_drude(a1)) {
10101 sprintf(msg,
"%s particle %d is bonded to non-parent atom %d",
10102 (is_lp(a2) ?
"Lone pair" :
"Drude"), a2+1, a1+1);
10105 atoms[a2].partner = a1;
10117 for(i=0; i<numAtoms; i++) {
10118 if ( ! hg[hg[i].GPID].isGP ) {
10120 sprintf(msg,
"child atom %d bonded only to child H atoms",i+1);
10123 if ( hg[i].isGP && is_hydrogen(i) ) {
10124 if ( hg[i].GPID == i )
continue;
10126 if ( is_hydrogen(hg[i].GPID) && hg[hg[i].GPID].GPID != i ) {
10128 sprintf(msg,
"H atom %d bonded only to child H atoms",i+1);
10134 if ( hg[i].atomsInGroup != 2 ) {
10136 sprintf(msg,
"H atom %d bonded to multiple H atoms",i+1);
10141 if ( hGPcount && ! CkMyPe() ) {
10142 iout <<
iWARN <<
"Found " << hGPcount <<
" H-H molecules.\n" <<
endi;
10146 for (i=0; i<numAtoms; ++i) {
10147 if ( hg[i].isGP ) hg[i].
GPID = i;
10155 for (i=0; i<numLphosts; ++i) {
10156 int a1 = lphosts[i].atom1;
10157 int a2 = lphosts[i].atom2;
10158 int a3 = lphosts[i].atom3;
10159 int a4 = lphosts[i].atom4;
10160 int m1 = hg[a1].
MPID;
10161 while ( hg[m1].MPID != m1 ) m1 = hg[m1].
MPID;
10162 int m2 = hg[a2].
MPID;
10163 while ( hg[m2].MPID != m2 ) m2 = hg[m2].
MPID;
10164 int m3 = hg[a3].
MPID;
10165 while ( hg[m3].MPID != m3 ) m3 = hg[m3].
MPID;
10166 int m4 = hg[a4].
MPID;
10167 while ( hg[m4].MPID != m4 ) m4 = hg[m4].
MPID;
10169 if ( m2 < mp ) mp = m2;
10170 if ( m3 < mp ) mp = m3;
10171 if ( m4 < mp ) mp = m4;
10179 for (i=0; i<numAtoms; ++i) {
10180 int mp = hg[i].
MPID;
10181 if ( hg[mp].MPID != mp ) {
10186 if ( allok )
break;
10188 for (i=0; i<numAtoms; ++i) {
10192 for (i=0; i<numAtoms; ++i) {
10198 for (i=0; i<numAtoms; i++) {
10209 numHydrogenGroups = 0;
10210 maxHydrogenGroupSize = 0;
10211 numMigrationGroups = 0;
10212 maxMigrationGroupSize = 0;
10213 for(i=0; i<numAtoms; i++)
10216 ++numMigrationGroups;
10218 if ( mgs > maxMigrationGroupSize ) maxMigrationGroupSize = mgs;
10221 ++numHydrogenGroups;
10223 if ( hgs > maxHydrogenGroupSize ) maxHydrogenGroupSize = hgs;
10227 hydrogenGroup.sort();
10232 for(i=0; i<numAtoms; ++i, --hgs) {
10234 if ( hg[i].isGP ) {
10236 parentid = hg[i].
atomID;
10239 sprintf(buff,
"Atom %d has bad hydrogen group size. " 10240 "Check for duplicate bonds.", parentid+1);
10244 if ( hg[i].isGP ) {
10246 sprintf(buff,
"Atom %d has bad hydrogen group size. " 10247 "Check for duplicate bonds.", parentid+1);
10255 for(i=0; i<numAtoms; ++i, --mgs) {
10257 if ( hg[i].isMP ) {
10259 parentid = hg[i].
atomID;
10262 sprintf(buff,
"Atom %d has bad migration group size.", parentid+1);
10266 if ( hg[i].isMP ) {
10268 sprintf(buff,
"Atom %d has bad migration group size.", parentid+1);
10276 for(i=0; i<numAtoms; i++) {
10277 atoms[hydrogenGroup[i].atomID].hydrogenList = i;
10283 for (i = 0; i < numAtoms; i++) {
10284 if (is_water(hg[i].atomID) && hg[i].isGP) {
10286 || ! is_drude(hg[i+1].atomID)
10287 || ! is_lp(hg[i+2].atomID)
10288 || ! is_hydrogen(hg[i+3].atomID)
10289 || ! is_hydrogen(hg[i+4].atomID) ) {
10291 sprintf(msg,
"Drude water molecule from HydrogenGroup i=%d " 10292 "starting at atom %d is not sorted\n", i, hg[i].atomID+1);
10299 else if (is_drude(hg[i].atomID)) {
10300 if (i < 1 || hg[i-1].atomID != hg[i].GPID) {
10302 sprintf(msg,
"Drude particle from HydrogenGroup i=%d must " 10303 "immediately follow its parent atom %d\n", i, hg[i].GPID+1);
10308 else if (is_lp(hg[i].atomID)) {
10310 sprintf(msg,
"Drude lonepair from HydrogenGroup i=%d " 10311 "at particle %d is NOT from water - unsupported\n",
10312 i, hg[i].atomID+1);
10322 for(i=0; i<numAtoms; i++)
10323 iout << i <<
" atomID=" << hydrogenGroup[i].atomID
10324 <<
" isGP=" << hydrogenGroup[i].isGP
10325 <<
" parent=" << hydrogenGroup[i].GPID
10326 <<
" #" << hydrogenGroup[i].atomsInGroup
10327 <<
" waterVal=" << hydrogenGroup[i].waterVal
10328 <<
" partner=" << atoms[i].partner
10329 <<
" hydrogenList=" << atoms[i].hydrogenList
10340 delete [] rigidBondLengths;
10341 rigidBondLengths =
new Real[numAtoms];
10342 if ( ! rigidBondLengths ) {
10343 NAMD_die(
"Memory allocation failed in Molecule::build_atom_status()\n");
10345 for (i=0; i<numAtoms; ++i) rigidBondLengths[i] = 0;
10350 for (i=0; i < numRealBonds; i++) {
10351 a1 = bonds[i].atom1;
10352 a2 = bonds[i].atom2;
10355 if (is_hydrogen(a2)) {
int tmp = a1; a1 = a2; a2 = tmp; }
10356 if (is_hydrogen(a1)) {
10357 if ( is_hydrogen(a2) ) {
10358 if ( ! is_water(a2) ) {
10359 rigidBondLengths[a1] = ( mode ==
RIGID_ALL ? x0 : 0. );
10360 rigidBondLengths[a2] = ( mode ==
RIGID_ALL ? x0 : 0. );
10362 }
else if ( is_water(a2) || mode ==
RIGID_ALL ) {
10363 rigidBondLengths[a1] = x0;
10364 if (is_water(a2)) r_oh = rigidBondLengths[a1];
10366 rigidBondLengths[a1] = 0.;
10371 if (is_lp(a2)) {
int tmp = a1; a1 = a2; a2 = tmp; }
10373 if (! is_water(a2) ) {
10377 sprintf(err_msg,
"ILLEGAL LONE PAIR AT INDEX %i\n" 10378 "LONE PAIRS ARE CURRENTLY ALLOWED ONLY ON WATER MOLECULES\n",
10382 rigidBondLengths[a1] = x0;
10391 int tmp = a1; a1 = a2; a2 = tmp;
10394 if (is_water(a2)) {
10401 sprintf(msg,
"ILLEGAL LONE PAIR AT INDEX %d\n" 10402 "LONE PAIRS ARE CURRENTLY ALLOWED ONLY ON WATER MOLECULES\n",
10412 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
10413 for( ; h_i != h_e; ++h_i ) {
10414 if ( h_i->
isGP ) rigidBondLengths[h_i->
atomID] = 0.;
10418 for (i=0; i < numAngles; i++) {
10419 a2 = angles[i].atom2;
10420 if ( ! is_water(a2) )
continue;
10421 if ( ! is_oxygen(a2) )
continue;
10422 a1 = angles[i].atom1;
10423 if ( ! is_hydrogen(a1) )
continue;
10424 a3 = angles[i].atom3;
10425 if ( ! is_hydrogen(a3) )
continue;
10426 if (is_lp(a2) || is_lp(a1) || is_lp(a3) ||
10427 is_drude(a2) || is_drude(a1) || is_drude(a3))
continue;
10428 if ( rigidBondLengths[a1] != rigidBondLengths[a3] ) {
10429 if (rigidBondLengths[a1] >0.3 && rigidBondLengths[a3] >0.3) {
10430 printf(
"length1: %f length2: %f\n", rigidBondLengths[a1], rigidBondLengths[a3]);
10432 NAMD_die(
"Asymmetric water molecule found??? This can't be right.\n");
10437 rigidBondLengths[a2] = 2. * rigidBondLengths[a1] * sin(0.5*t0);
10438 r_hh = rigidBondLengths[a2];
10442 int numBondWaters = 0;
10443 int numFailedWaters = 0;
10445 for (i=0; i < numRealBonds; i++) {
10446 a1 = bonds[i].atom1;
10447 a2 = bonds[i].atom2;
10448 if ( ! is_hydrogen(a1) )
continue;
10449 if ( ! is_hydrogen(a2) )
continue;
10450 int ma1 = get_mother_atom(a1);
10451 int ma2 = get_mother_atom(a2);
10452 if ( ma1 != ma2 )
continue;
10453 if ( ! is_water(ma1) )
continue;
10454 if ( rigidBondLengths[ma1] != 0. )
continue;
10457 rigidBondLengths[ma1] = x0;
10464 if (r_oh < 0.0 || r_hh < 0.0) {
10466 NAMD_die(
"Failed to find water bond lengths\n");
10468 r_ohc = sqrt(r_oh * r_oh - 0.25 * r_hh * r_hh);
10472 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
10473 for( ; h_i != h_e; ++h_i ) {
10475 rigidBondLengths[h_i->
atomID] == 0. ) {
10476 if ( h_i + 1 == h_e || h_i + 2 == h_e ||
10477 h_i[1].isGP || h_i[2].isGP || h_i->
atomsInGroup != 3 ) {
10478 NAMD_die(
"Abnormal water detected.");
10480 if ( CkNumNodes() > 1 ) {
10481 NAMD_die(
"Unable to determine H-H distance for rigid water because structure has neither H-O-H angle nor H-H bond.");
10487 get_atomtype(btmp.
atom1),
10488 get_atomtype(btmp.
atom2),
10494 rigidBondLengths[h_i->
atomID] = x0;
10501 if ( numBondWaters + numFailedWaters ) {
10502 iout <<
iWARN <<
"Missing angles for " <<
10503 ( numBondWaters + numFailedWaters ) <<
" waters.\n" <<
endi;
10505 if ( numBondWaters ) {
10506 iout <<
iWARN <<
"Obtained H-H distance from bond parameters for " <<
10507 numBondWaters <<
" waters.\n" <<
endi;
10508 iout <<
iWARN <<
"This would not be possible in a multi-process run.\n" <<
endi;
10510 if ( numFailedWaters ) {
10511 iout <<
iERROR <<
"Failed to obtain H-H distance from angles or bonds for " <<
10512 numFailedWaters <<
" waters.\n" <<
endi;
10520 for (i=0; i<numAtoms; ++i) rigidBondLengths[i] *= -1;
10522 for (i=0; i<numAtoms; ++i) {
10523 if ( ! is_water(i) ) rigidBondLengths[i] *= -1;
10529 for (i=0; i<numAtoms; ++i) {
10530 if ( rigidBondLengths[i] > 0. ) ++numRigidBonds;
10551 BigReal LJAvgA, LJAvgB, LJAvgA1, LJAvgB1, LJAvgA2, LJAvgB2;
10552 int numLJsites, numLJsites1, numLJsites2;
10564 const int ss_dim = soluteScalingOn ? ss_num_vdw_params : 0;
10572 if (soluteScalingOn && (ss_vdw_type == NULL)) {
10573 NAMD_bug(
"Solute scaling is used but ss_vdw_type is NULL. " 10574 "It is likely that compute_LJcorrection() is called " 10575 "before build_ss_flags().");
10577 Real A, B, A14, B14;
10578 Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
10579 Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
10580 Real *ATable =
new Real[LJtypecount*LJtypecount];
10581 Real *BTable =
new Real[LJtypecount*LJtypecount];
10582 int useGeom =
simParams->vdwGeometricSigma;
10584 for (
int i = 0; i < LJtypecount; i++) {
10585 for (
int j = 0; j < LJtypecount; j++) {
10586 const int i_type = soluteScalingOn ? ((i >= table_dim_org) ? ss_vdw_type[i-table_dim_org] : i) : i;
10587 const int j_type = soluteScalingOn ? ((j >= table_dim_org) ? ss_vdw_type[j-table_dim_org] : j) : j;
10589 ATable[i*LJtypecount + j] = A;
10590 BTable[i*LJtypecount + j] = B;
10593 params->
get_vdw_params(&sigma_i,&epsilon_i,&sigma_i14,&epsilon_i14,i_type);
10594 params->
get_vdw_params(&sigma_j,&epsilon_j,&sigma_j14,&epsilon_j14,j_type);
10596 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10597 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10598 sigma_ij *= sigma_ij*sigma_ij;
10599 sigma_ij *= sigma_ij;
10600 ATable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10601 BTable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij;
10606 int *numAtomsByLjType =
new int[LJtypecount];
10607 for (
int i=0; i < LJtypecount; i++) {numAtomsByLjType[i]=0;}
10608 for (
int i=0; i < numAtoms; i++) {numAtomsByLjType[atoms[i].vdw_type]++;}
10615 for (
int i=0; i < LJtypecount; i++) {
10616 for (
int j=0; j < LJtypecount; j++) {
10617 A = ATable[i*LJtypecount + j];
10618 B = BTable[i*LJtypecount + j];
10619 if (!A && !B)
continue;
10620 npairs = (numAtomsByLjType[i] - int(i==j))*
BigReal(numAtomsByLjType[j]);
10621 sumOfAs += npairs*A;
10622 sumOfBs += npairs*B;
10624 if (i==j) numLJsites += numAtomsByLjType[i];
10627 delete [] numAtomsByLjType;
10631 LJAvgA = sumOfAs / count;
10632 LJAvgB = sumOfBs / count;
10653 numLJsites1 = numLJsites2 = numLJsites;
10654 int alch_counter = 0;
10655 for (
int i=0; i < numAtoms; ++i) {
10657 if (get_fep_type(i) == 2 || get_fep_type(i) == 4) alchFlagi = -1;
10658 if (get_fep_type(i) == 1 || get_fep_type(i) == 3) alchFlagi = 1;
10660 &A, &B, &A14, &B14)) {
10664 &epsilon_i14, atoms[i].vdw_type);
10666 useGeom ? sqrt(sigma_i*sigma_i) : 0.5*(sigma_i+sigma_i);
10667 BigReal epsilon_ii = sqrt(epsilon_i*epsilon_i);
10669 sigma_ii *= sigma_ii*sigma_ii;
10670 sigma_ii *= sigma_ii;
10671 A = 4.0*sigma_ii*epsilon_ii*sigma_ii;
10672 B = 4.0*sigma_ii*epsilon_ii;
10675 if (alchFlagi == 1) numLJsites2--;
10676 else if (alchFlagi == -1) numLJsites1--;
10678 for (
int j=i+1; j < numAtoms; ++j) {
10680 if (get_fep_type(j) == 2 || get_fep_type(j) == 4) alchFlagj = -1;
10681 if (get_fep_type(j) == 1 || get_fep_type(j) == 3) alchFlagj = 1;
10682 int alchFlagSum = alchFlagi + alchFlagj;
10685 if (alchFlagi == 0 && alchFlagj == 0)
continue;
10688 &A, &B, &A14, &B14)) {
10692 &epsilon_i14, atoms[i].vdw_type);
10694 &epsilon_j14, atoms[j].vdw_type);
10696 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10697 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10699 sigma_ij *= sigma_ij*sigma_ij;
10700 sigma_ij *= sigma_ij;
10701 A = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10702 B = 4.0*sigma_ij*epsilon_ij;
10704 if (!A && !B)
continue;
10709 if ( alchFlagSum > 0 ){
10714 else if ( alchFlagSum < 0 ){
10730 if ( alchFlagi == 1 || alchFlagi == -1 ) alch_counter++;
10731 if ( alch_counter == (numFepInitial + numFepFinal) )
break;
10733 LJAvgA = sumOfAs / count;
10734 LJAvgB = sumOfBs / count;
10736 LJAvgA1 = sumOfAs1 / count1;
10737 LJAvgB1 = sumOfBs1 / count1;
10743 LJAvgA2 = sumOfAs2 / count2;
10744 LJAvgB2 = sumOfBs2 / count2;
10749 if ( ! CkMyPe() ) {
10750 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 10751 <<
"ENERGY AND PRESSURE\n" <<
endi;
10752 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A0 AND B0 COEFFICIENTS " 10753 << LJAvgA2 <<
" AND " << LJAvgB2 <<
"\n" <<
endi;
10754 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A1 AND B1 COEFFICIENTS " 10755 << LJAvgA1 <<
" AND " << LJAvgB1 <<
"\n" <<
endi;
10757 numLJsites = (numLJsites1 + numLJsites2 - numLJsites);
10758 LJAvgA1 *=
BigReal(numLJsites1)*numLJsites1;
10759 LJAvgB1 *=
BigReal(numLJsites1)*numLJsites1;
10760 LJAvgA2 *=
BigReal(numLJsites2)*numLJsites2;
10761 LJAvgB2 *=
BigReal(numLJsites2)*numLJsites2;
10764 LJAvgA1 = LJAvgB1 = LJAvgA2 = LJAvgB2 = 0;
10766 if ( ! CkMyPe() ) {
10767 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 10768 <<
"ENERGY AND PRESSURE\n" <<
endi;
10769 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A AND B COEFFICIENTS " 10770 << LJAvgA <<
" AND " << LJAvgB <<
"\n" <<
endi;
10773 LJAvgA *=
BigReal(numLJsites)*numLJsites;
10774 LJAvgB *=
BigReal(numLJsites)*numLJsites;
10783 BigReal rswitch2 = rswitch*rswitch;
10784 BigReal rswitch3 = rswitch*rswitch2;
10785 BigReal rswitch4 = rswitch2*rswitch2;
10786 BigReal rswitch5 = rswitch2*rswitch3;
10787 BigReal rswitch6 = rswitch3*rswitch3;
10813 BigReal int_U_gofr_A, int_rF_gofr_A, int_U_gofr_B, int_rF_gofr_B;
10832 BigReal rsum3 = (rcut + rswitch)*(rcut + rswitch)*(rcut + rswitch);
10833 int_U_gofr_A = int_rF_gofr_A = (16*
PI*(3*rcut4 + 9*rcut3*rswitch
10834 + 11*rcut2*rswitch2 + 9*rcut*rswitch3
10836 / (315*rcut5*rswitch5*rsum3));
10837 int_U_gofr_B = int_rF_gofr_B = -16*
PI / (3*rsum3);
10864 BigReal lnr = log(rcut/rswitch);
10888 int_U_gofr_A = int_rF_gofr_A =\
10889 16*
PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
10890 int_U_gofr_B = int_rF_gofr_B = -4*
PI*lnr / (rcut3 - rswitch3);
10927 int_rF_gofr_A = 8*
PI / (9*rcut9);
10928 int_rF_gofr_B = -4*
PI / (3*rcut3);
10929 int_U_gofr_A = 2*
PI / (9*rcut9);
10930 int_U_gofr_B = -2*
PI / (3*rcut3);
10933 tail_corr_virial = int_rF_gofr_A*LJAvgA + int_rF_gofr_B*LJAvgB;
10934 tail_corr_ener = int_U_gofr_A*LJAvgA + int_U_gofr_B*LJAvgB;
10936 tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1 -
10938 tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1 -
10940 tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2 -
10942 tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2 -
10961 BigReal LJAvgA, LJAvgB, LJAvgA1, LJAvgB1, LJAvgA2, LJAvgB2;
10962 int numLJsites, numLJsites1, numLJsites2;
10974 const int ss_dim = soluteScalingOn ? ss_num_vdw_params : 0;
10982 if (soluteScalingOn && (ss_vdw_type == NULL)) {
10983 NAMD_bug(
"Solute scaling is used but ss_vdw_type is NULL. " 10984 "It is likely that compute_LJcorrection_alternative() is called " 10985 "before build_ss_flags().");
10987 Real A, B, A14, B14;
10988 Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
10989 Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
10990 Real *ATable =
new Real[LJtypecount*LJtypecount];
10991 Real *BTable =
new Real[LJtypecount*LJtypecount];
10992 int useGeom =
simParams->vdwGeometricSigma;
10994 for (
int i = 0; i < LJtypecount; i++) {
10995 for (
int j = 0; j < LJtypecount; j++) {
10996 const int i_type = soluteScalingOn ? ((i >= table_dim_org) ? ss_vdw_type[i-table_dim_org] : i) : i;
10997 const int j_type = soluteScalingOn ? ((j >= table_dim_org) ? ss_vdw_type[j-table_dim_org] : j) : j;
10999 ATable[i*LJtypecount + j] = A;
11000 BTable[i*LJtypecount + j] = B;
11003 params->
get_vdw_params(&sigma_i,&epsilon_i,&sigma_i14,&epsilon_i14,i_type);
11004 params->
get_vdw_params(&sigma_j,&epsilon_j,&sigma_j14,&epsilon_j14,j_type);
11006 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
11007 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
11008 sigma_ij *= sigma_ij*sigma_ij;
11009 sigma_ij *= sigma_ij;
11010 ATable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij*sigma_ij;
11011 BTable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij;
11016 int *numAtomsByLjType =
new int[LJtypecount];
11017 for (
int i=0; i < LJtypecount; i++) {numAtomsByLjType[i]=0;}
11018 for (
int i=0; i < numAtoms; i++) {numAtomsByLjType[atoms[i].vdw_type]++;}
11023 for (
int i=0; i < LJtypecount; i++) {
11024 for (
int j=0; j < LJtypecount; j++) {
11025 A = ATable[i*LJtypecount + j];
11026 B = BTable[i*LJtypecount + j];
11027 npairs =
BigReal(numAtomsByLjType[i])*
BigReal(numAtomsByLjType[j]);
11028 sumOfAs += npairs*A;
11029 sumOfBs += npairs*B;
11032 delete [] numAtomsByLjType;
11056 int alch_counter = 0;
11057 for (
int i=0; i < numAtoms; ++i) {
11059 if (get_fep_type(i) == 2 || get_fep_type(i) == 4) alchFlagi = -1;
11060 if (get_fep_type(i) == 1 || get_fep_type(i) == 3) alchFlagi = 1;
11062 for (
int j=i; j < numAtoms; ++j) {
11064 if (get_fep_type(j) == 2 || get_fep_type(j) == 4) alchFlagj = -1;
11065 if (get_fep_type(j) == 1 || get_fep_type(j) == 3) alchFlagj = 1;
11066 int alchFlagSum = alchFlagi + alchFlagj;
11069 if (alchFlagi == 0 && alchFlagj == 0)
continue;
11072 &A, &B, &A14, &B14)) {
11076 &epsilon_i14, atoms[i].vdw_type);
11078 &epsilon_j14, atoms[j].vdw_type);
11080 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
11081 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
11083 sigma_ij *= sigma_ij*sigma_ij;
11084 sigma_ij *= sigma_ij;
11085 A = 4.0*sigma_ij*epsilon_ij*sigma_ij;
11086 B = 4.0*sigma_ij*epsilon_ij;
11091 if ( alchFlagSum > 0 ){
11094 }
else if ( alchFlagSum < 0 ){
11106 if ( alchFlagi == 1 || alchFlagi == -1 ) alch_counter++;
11107 if ( alch_counter == (numFepInitial + numFepFinal) )
break;
11112 LJAvgA1 = sumOfAs1;
11113 LJAvgB1 = sumOfBs1;
11114 LJAvgA2 = sumOfAs2;
11115 LJAvgB2 = sumOfBs2;
11117 if ( ! CkMyPe() ) {
11118 BigReal inv_sizeSq = 1.0 / (numAtoms * numAtoms);
11119 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 11120 <<
"ENERGY AND PRESSURE\n" <<
endi;
11121 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A0 AND B0 COEFFICIENTS " 11122 << LJAvgA2*inv_sizeSq <<
" AND " << LJAvgB2*inv_sizeSq <<
"\n" <<
endi;
11123 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A1 AND B1 COEFFICIENTS " 11124 << LJAvgA1*inv_sizeSq <<
" AND " << LJAvgB1*inv_sizeSq <<
"\n" <<
endi;
11127 LJAvgA1 = LJAvgB1 = LJAvgA2 = LJAvgB2 = 0;
11129 if ( ! CkMyPe() ) {
11130 BigReal inv_sizeSq = 1.0 / (numAtoms * numAtoms);
11131 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 11132 <<
"ENERGY AND PRESSURE\n" <<
endi;
11133 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A AND B COEFFICIENTS " 11134 << LJAvgA*inv_sizeSq <<
" AND " << LJAvgB*inv_sizeSq <<
"\n" <<
endi;
11145 BigReal rswitch2 = rswitch*rswitch;
11146 BigReal rswitch3 = rswitch*rswitch2;
11147 BigReal rswitch4 = rswitch2*rswitch2;
11148 BigReal rswitch5 = rswitch2*rswitch3;
11149 BigReal rswitch6 = rswitch3*rswitch3;
11175 BigReal int_U_gofr_A, int_rF_gofr_A, int_U_gofr_B, int_rF_gofr_B;
11194 BigReal rsum3 = (rcut + rswitch)*(rcut + rswitch)*(rcut + rswitch);
11195 int_U_gofr_A = int_rF_gofr_A = (16*
PI*(3*rcut4 + 9*rcut3*rswitch
11196 + 11*rcut2*rswitch2 + 9*rcut*rswitch3
11198 / (315*rcut5*rswitch5*rsum3));
11199 int_U_gofr_B = int_rF_gofr_B = -16*
PI / (3*rsum3);
11226 BigReal lnr = log(rcut/rswitch);
11250 int_U_gofr_A = int_rF_gofr_A =\
11251 16*
PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
11252 int_U_gofr_B = int_rF_gofr_B = -4*
PI*lnr / (rcut3 - rswitch3);
11289 int_rF_gofr_A = 8*
PI / (9*rcut9);
11290 int_rF_gofr_B = -4*
PI / (3*rcut3);
11291 int_U_gofr_A = 2*
PI / (9*rcut9);
11292 int_U_gofr_B = -2*
PI / (3*rcut3);
11295 tail_corr_virial = int_rF_gofr_A*LJAvgA + int_rF_gofr_B*LJAvgB;
11296 tail_corr_ener = int_U_gofr_A*LJAvgA + int_U_gofr_B*LJAvgB;
11298 tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1 -
11300 tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1 -
11302 tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2 -
11304 tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2 -
11316 const BigReal corr = (doti ? 0.0 : tail_corr_ener);
11317 return scl*(corr + vdw_lambda_1*tail_corr_dUdl_1 +
11318 vdw_lambda_2*tail_corr_dUdl_2);
11321 return scl*tail_corr_ener;
11331 return scl*(tail_corr_virial + vdw_lambda_1*tail_corr_virial_1 +
11332 vdw_lambda_2*tail_corr_virial_2);
11335 return scl*tail_corr_virial;
11340 #ifdef MEM_OPT_VERSION 11343 int Molecule::checkExclByIdx(
int idx1,
int atom1,
int atom2)
const {
11345 int amin = exclChkSigPool[idx1].min;
11346 int amax = exclChkSigPool[idx1].max;
11347 int dist21 = atom2 - atom1;
11348 if ( dist21 < amin || dist21 > amax )
return 0;
11349 else return exclChkSigPool[idx1].flags[dist21-amin];
11355 int amin = all_exclusions[atom1].min;
11356 int amax = all_exclusions[atom1].max;
11357 if ( atom2 < amin || atom2 > amax )
return 0;
11358 else return all_exclusions[atom1].flags[atom2-amin];
11378 is_lonepairs_psf = 0;
11386 read_parm(amber_data);
11388 #ifndef MEM_OPT_VERSION 11391 assignLCPOTypes( 1 );
11403 is_lonepairs_psf = 0;
11411 read_parm(amber_data);
11413 #ifndef MEM_OPT_VERSION 11416 assignLCPOTypes( 1 );
11421 #ifdef MEM_OPT_VERSION 11422 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
11424 int i,j,ntheth,nphih,current_index,a1,a2,
11425 max,min,index,found;
11428 NAMD_die(
"No data read from parm file yet!");
11431 numAtoms = amber_data->
Natom;
11432 atoms =
new Atom[numAtoms];
11444 if (atoms == NULL || atomNames == NULL )
11445 NAMD_die(
"memory allocation failed when reading atom information");
11448 for (i = 0; i < numAtoms; ++i) {
11449 const int resname_length = amber_data->
ResNames[amber_data->
AtomRes[i]].size();
11450 const int atomname_length = amber_data->
AtomNames[i].size();
11451 const int atomtype_length = amber_data->
AtomSym[i].size();
11452 atomNames[i].resname = nameArena->getNewArray(resname_length+1);
11453 atomNames[i].atomname = nameArena->getNewArray(atomname_length+1);
11454 atomNames[i].atomtype = nameArena->getNewArray(atomtype_length+1);
11455 if (atomNames[i].resname == NULL || atomNames[i].atomname == NULL || atomNames[i].atomtype == NULL)
11456 NAMD_die(
"memory allocation failed when reading atom information");
11458 strncpy(atomNames[i].resname, amber_data->
ResNames[amber_data->
AtomRes[i]].c_str(), resname_length+1);
11459 strncpy(atomNames[i].atomname, amber_data->
AtomNames[i].c_str(), atomname_length+1);
11460 strncpy(atomNames[i].atomtype, amber_data->
AtomSym[i].c_str(), atomtype_length+1);
11462 strtok(atomNames[i].resname,
" ");
11463 strtok(atomNames[i].atomname,
" ");
11464 strtok(atomNames[i].atomtype,
" ");
11465 atoms[i].mass = amber_data->
Masses[i];
11512 atoms[i].charge = amber_data->
Charges[i] / 18.2223;
11513 atoms[i].vdw_type = amber_data->
Iac[i] - 1;
11515 if ( tmpResLookup ) tmpResLookup =
11518 if(atomSegResids) {
11520 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
11527 }
else if (atoms[i].mass <= 0.05) {
11529 ++numZeroMassAtoms;
11530 }
else if (atoms[i].mass < 1.0) {
11536 }
else if (atoms[i].mass <=3.5) {
11538 }
else if ((atomNames[i].atomname[0] ==
'O') &&
11539 (atoms[i].mass >= 14.0) &&
11540 (atoms[i].mass <= 18.0)) {
11553 if (amber_data->
Nbonh + amber_data->
Nbona > 0) {
11555 if (bonds == NULL || amber_data->
Nbona < 0)
11556 NAMD_die(
"memory allocation failed when reading bond information");
11558 for (i=0; i<amber_data->
Nbonh; ++i)
11559 { bonds[numBonds].atom1 = amber_data->
BondHAt1[i] / 3;
11560 bonds[numBonds].atom2 = amber_data->
BondHAt2[i] / 3;
11561 bonds[numBonds].bond_type = amber_data->
BondHNum[i] - 1;
11562 if (bonds[numBonds].atom1>=numAtoms || bonds[numBonds].atom2>=numAtoms ||
11563 bonds[numBonds].bond_type>=amber_data->
Numbnd)
11564 {
char err_msg[128];
11565 sprintf(err_msg,
"BOND (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11573 for (i=amber_data->
Nbonh; i<amber_data->Nbonh+amber_data->
Nbona; ++i)
11574 { bonds[numBonds].atom1 = amber_data->
BondAt1[i-amber_data->
Nbonh] / 3;
11575 bonds[numBonds].atom2 = amber_data->
BondAt2[i-amber_data->
Nbonh] / 3;
11576 bonds[numBonds].bond_type = amber_data->
BondNum[i-amber_data->
Nbonh] - 1;
11577 if (bonds[i].atom1>=numAtoms || bonds[i].atom2>=numAtoms ||
11578 bonds[i].bond_type>=amber_data->
Numbnd)
11579 {
char err_msg[128];
11580 sprintf(err_msg,
"BOND (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-amber_data->
Nbonh);
11589 if ( numBonds != amber_data->
Nbonh + amber_data->
Nbona) {
11591 " bonds with zero force constants.\n" <<
endi;
11593 "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
endi;
11599 { ntheth = amber_data->
Ntheth;
11600 angles =
new Angle[numAngles];
11601 if (angles == NULL || numAngles < ntheth)
11602 NAMD_die(
"memory allocation failed when reading angle information");
11604 for (i=0; i<ntheth; ++i)
11605 { angles[i].atom1 = amber_data->
AngleHAt1[i] / 3;
11606 angles[i].atom2 = amber_data->
AngleHAt2[i] / 3;
11607 angles[i].atom3 = amber_data->
AngleHAt3[i] / 3;
11608 angles[i].angle_type = amber_data->
AngleHNum[i] - 1;
11609 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11610 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11611 {
char err_msg[128];
11612 sprintf(err_msg,
"ANGLE (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11617 for (i=ntheth; i<numAngles; ++i)
11618 { angles[i].atom1 = amber_data->
AngleAt1[i-ntheth] / 3;
11619 angles[i].atom2 = amber_data->
AngleAt2[i-ntheth] / 3;
11620 angles[i].atom3 = amber_data->
AngleAt3[i-ntheth] / 3;
11621 angles[i].angle_type = amber_data->
AngleNum[i-ntheth] - 1;
11622 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11623 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11624 {
char err_msg[128];
11625 sprintf(err_msg,
"ANGLE (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-ntheth);
11640 if (exclusions == NULL && amber_data->
Nnb > 0)
11641 NAMD_die(
"memory allocation failed when reading exclusion information");
11643 for (i=0; i<numAtoms; ++i)
11644 for (j=0; j<amber_data->
Iblo[i]; ++j)
11645 {
if (current_index >= amber_data->
Nnb)
11646 {
char err_msg[128];
11647 sprintf(err_msg,
"EXCLUSION INDEX EXCEEDS NUMBER OF EXLCUSIONS %d IN AMBER FILE, AT ATOM #%d\n",
11648 amber_data->
Nnb, i+1);
11653 if (amber_data->
ExclAt[current_index] != 0)
11656 a2 = amber_data->
ExclAt[current_index] - 1;
11662 sprintf(err_msg,
"Atom #%d has exclusion with atom #%d, in reverse order.", i+1, a2+1);
11666 {
char err_msg[128];
11667 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN AMBER FILE\n", i+1);
11670 else if (a2 >= numAtoms)
11671 {
char err_msg[128];
11672 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN AMBER FILE",
11673 a2+1, numAtoms, current_index+1);
11676 exclusions[numExclusions].atom1 = i;
11677 exclusions[numExclusions].atom2 = a2;
11682 if (current_index < amber_data->Nnb)
11683 {
char err_msg[128];
11684 sprintf(err_msg,
"Num of exclusions recorded (%d) is smaller than what it's supposed to be (%d)",
11685 current_index,amber_data->
Nnb);
11691 numDihedrals = amber_data->
Nphih + amber_data->
Nphia;
11692 if (numDihedrals > 0)
11693 { nphih = amber_data->
Nphih;
11694 dihedrals =
new Dihedral[numDihedrals];
11695 if (dihedrals == NULL || numDihedrals < nphih)
11696 NAMD_die(
"memory allocation failed when reading dihedral information");
11698 for (i=0; i<nphih; ++i)
11699 { dihedrals[i].atom1 = amber_data->
DihHAt1[i] / 3;
11700 dihedrals[i].atom2 = amber_data->
DihHAt2[i] / 3;
11701 dihedrals[i].atom3 = amber_data->
DihHAt3[i] / 3;
11702 dihedrals[i].atom4 = amber_data->
DihHAt4[i] / 3;
11703 dihedrals[i].dihedral_type = amber_data->
DihHNum[i] - 1;
11706 for (i=nphih; i<numDihedrals; ++i)
11707 { dihedrals[i].atom1 = amber_data->
DihAt1[i-nphih] / 3;
11708 dihedrals[i].atom2 = amber_data->
DihAt2[i-nphih] / 3;
11709 dihedrals[i].atom3 = amber_data->
DihAt3[i-nphih] / 3;
11710 dihedrals[i].atom4 = amber_data->
DihAt4[i-nphih] / 3;
11711 dihedrals[i].dihedral_type = amber_data->
DihNum[i-nphih] - 1;
11726 for (i=0; i<numDihedrals; ++i)
11727 {
if (dihedrals[i].atom3 < 0 || dihedrals[i].atom4 < 0)
11728 { dihedrals[i].atom3 = abs(dihedrals[i].atom3);
11729 dihedrals[i].atom4 = abs(dihedrals[i].atom4);
11732 {
if (dihedrals[i].atom1 < dihedrals[i].atom4)
11733 a1=dihedrals[i].atom1, a2=dihedrals[i].atom4;
11735 a1=dihedrals[i].atom4, a2=dihedrals[i].atom1;
11739 min=0, max=numExclusions-1;
11740 while (!found && min<=max)
11741 { index = (min+max)/2;
11742 if (exclusions[index].atom1 == a1)
11744 else if (exclusions[index].atom1 < a1)
11750 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11753 for (j=index-1; j>=0 && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; --j);
11754 if (j<0 || exclusions[j].atom1!=a1)
11755 for (j=index; j<numExclusions && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; ++j);
11756 if (j<numExclusions && exclusions[j].atom1==a1)
11757 exclusions[j].modified = 1;
11759 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11761 if (dihedrals[i].atom1>=numAtoms || dihedrals[i].atom2>=numAtoms ||
11762 dihedrals[i].atom3>=numAtoms || dihedrals[i].atom4>=numAtoms ||
11763 dihedrals[i].dihedral_type>=amber_data->
Nptra)
11764 {
char err_msg[128];
11765 sprintf(err_msg,
"DIHEDRAL # %d OVERFLOW IN PARM FILE", i+1);
11771 if (amber_data -> HasCMAP) {
11773 if (numCrossterms > 0) {
11774 crossterms =
new Crossterm[numCrossterms];
11776 for (i = 0; i < numCrossterms; ++i) {
11779 crossterms[i].atom2 = amber_data->
CMAPIndex[6*i+1]-1;
11780 crossterms[i].atom3 = amber_data->
CMAPIndex[6*i+2]-1;
11781 crossterms[i].atom4 = amber_data->
CMAPIndex[6*i+3]-1;
11783 crossterms[i].atom5 = amber_data->
CMAPIndex[6*i+1]-1;
11784 crossterms[i].atom6 = amber_data->
CMAPIndex[6*i+2]-1;
11785 crossterms[i].atom7 = amber_data->
CMAPIndex[6*i+3]-1;
11786 crossterms[i].atom8 = amber_data->
CMAPIndex[6*i+4]-1;
11788 crossterms[i].crossterm_type = amber_data->
CMAPIndex[6*i+5]-1;
11793 numRealBonds = numBonds;
11794 build_atom_status();
11810 void Molecule::read_parm(
Ambertoppar *amber_data)
11812 #ifdef MEM_OPT_VERSION 11813 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
11815 int i,j,ntheth,nphih,current_index,a1,a2,
11816 max,min,index,found;
11819 NAMD_die(
"No data read from parm file yet!");
11822 numAtoms = amber_data->
Natom;
11823 atoms =
new Atom[numAtoms];
11830 if (atoms == NULL || atomNames == NULL )
11831 NAMD_die(
"memory allocation failed when reading atom information");
11833 for (i=0; i<numAtoms; ++i)
11834 { atomNames[i].resname = nameArena->getNewArray(5);
11835 atomNames[i].atomname = nameArena->getNewArray(5);
11836 atomNames[i].atomtype = nameArena->getNewArray(5);
11837 if (atomNames[i].resname == NULL || atomNames[i].atomname == NULL || atomNames[i].atomtype == NULL)
11838 NAMD_die(
"memory allocation failed when reading atom information");
11839 for (j=0; j<4; ++j)
11840 { atomNames[i].resname[j] = amber_data->
ResNames[amber_data->
AtomRes[i]*4+j];
11841 atomNames[i].atomname[j] = amber_data->
AtomNames[i*4+j];
11842 atomNames[i].atomtype[j] = amber_data->
AtomSym[i*4+j];
11844 atomNames[i].resname[4] = atomNames[i].atomname[4] = atomNames[i].atomtype[4] =
'\0';
11845 strtok(atomNames[i].resname,
" ");
11846 strtok(atomNames[i].atomname,
" ");
11847 strtok(atomNames[i].atomtype,
" ");
11848 atoms[i].mass = amber_data->
Masses[i];
11850 atoms[i].charge = amber_data->
Charges[i] / 18.2223;
11851 atoms[i].vdw_type = amber_data->
Iac[i] - 1;
11854 if ( tmpResLookup ) tmpResLookup =
11857 if(atomSegResids) {
11859 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
11867 }
else if (atoms[i].mass <= 0.05) {
11868 ++numZeroMassAtoms;
11870 }
else if (atoms[i].mass < 1.0) {
11872 }
else if (atoms[i].mass <=3.5) {
11874 }
else if ((atomNames[i].atomname[0] ==
'O') &&
11875 (atoms[i].mass >= 14.0) &&
11876 (atoms[i].mass <= 18.0)) {
11889 if (amber_data->
Nbonh + amber_data->
Nbona > 0)
11891 if (bonds == NULL || amber_data->
Nbona < 0)
11892 NAMD_die(
"memory allocation failed when reading bond information");
11894 for (i=0; i<amber_data->
Nbonh; ++i)
11895 { bonds[numBonds].atom1 = amber_data->
BondHAt1[i] / 3;
11896 bonds[numBonds].atom2 = amber_data->
BondHAt2[i] / 3;
11897 bonds[numBonds].bond_type = amber_data->
BondHNum[i] - 1;
11898 if (bonds[numBonds].atom1>=numAtoms || bonds[numBonds].atom2>=numAtoms ||
11899 bonds[numBonds].bond_type>=amber_data->
Numbnd)
11900 {
char err_msg[128];
11901 sprintf(err_msg,
"BOND (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11909 for (i=amber_data->
Nbonh; i<amber_data->Nbonh+amber_data->
Nbona; ++i)
11910 { bonds[numBonds].atom1 = amber_data->
BondAt1[i-amber_data->
Nbonh] / 3;
11911 bonds[numBonds].atom2 = amber_data->
BondAt2[i-amber_data->
Nbonh] / 3;
11912 bonds[numBonds].bond_type = amber_data->
BondNum[i-amber_data->
Nbonh] - 1;
11913 if (bonds[i].atom1>=numAtoms || bonds[i].atom2>=numAtoms ||
11914 bonds[i].bond_type>=amber_data->
Numbnd)
11915 {
char err_msg[128];
11916 sprintf(err_msg,
"BOND (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-amber_data->
Nbonh);
11925 if ( numBonds != amber_data->
Nbonh + amber_data->
Nbona) {
11927 " bonds with zero force constants.\n" <<
endi;
11929 "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
endi;
11935 { ntheth = amber_data->
Ntheth;
11936 angles =
new Angle[numAngles];
11937 if (angles == NULL || numAngles < ntheth)
11938 NAMD_die(
"memory allocation failed when reading angle information");
11940 for (i=0; i<ntheth; ++i)
11941 { angles[i].atom1 = amber_data->
AngleHAt1[i] / 3;
11942 angles[i].atom2 = amber_data->
AngleHAt2[i] / 3;
11943 angles[i].atom3 = amber_data->
AngleHAt3[i] / 3;
11944 angles[i].angle_type = amber_data->
AngleHNum[i] - 1;
11945 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11946 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11947 {
char err_msg[128];
11948 sprintf(err_msg,
"ANGLE (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11953 for (i=ntheth; i<numAngles; ++i)
11954 { angles[i].atom1 = amber_data->
AngleAt1[i-ntheth] / 3;
11955 angles[i].atom2 = amber_data->
AngleAt2[i-ntheth] / 3;
11956 angles[i].atom3 = amber_data->
AngleAt3[i-ntheth] / 3;
11957 angles[i].angle_type = amber_data->
AngleNum[i-ntheth] - 1;
11958 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11959 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11960 {
char err_msg[128];
11961 sprintf(err_msg,
"ANGLE (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-ntheth);
11976 if (exclusions == NULL && amber_data->
Nnb > 0)
11977 NAMD_die(
"memory allocation failed when reading exclusion information");
11979 for (i=0; i<numAtoms; ++i)
11980 for (j=0; j<amber_data->
Iblo[i]; ++j)
11981 {
if (current_index >= amber_data->
Nnb)
11982 {
char err_msg[128];
11983 sprintf(err_msg,
"EXCLUSION INDEX EXCEEDS NUMBER OF EXLCUSIONS %d IN AMBER FILE, AT ATOM #%d\n",
11984 amber_data->
Nnb, i+1);
11989 if (amber_data->
ExclAt[current_index] != 0)
11992 a2 = amber_data->
ExclAt[current_index] - 1;
11998 sprintf(err_msg,
"Atom #%d has exclusion with atom #%d, in reverse order.", i+1, a2+1);
12002 {
char err_msg[128];
12003 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN AMBER FILE\n", i+1);
12006 else if (a2 >= numAtoms)
12007 {
char err_msg[128];
12008 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN AMBER FILE",
12009 a2+1, numAtoms, current_index+1);
12012 exclusions[numExclusions].atom1 = i;
12013 exclusions[numExclusions].atom2 = a2;
12018 if (current_index < amber_data->Nnb)
12019 {
char err_msg[128];
12020 sprintf(err_msg,
"Num of exclusions recorded (%d) is smaller than what it's supposed to be (%d)",
12021 current_index,amber_data->
Nnb);
12027 numDihedrals = amber_data->
Nphih + amber_data->
Nphia;
12028 if (numDihedrals > 0)
12029 { nphih = amber_data->
Nphih;
12030 dihedrals =
new Dihedral[numDihedrals];
12031 if (dihedrals == NULL || numDihedrals < nphih)
12032 NAMD_die(
"memory allocation failed when reading dihedral information");
12034 for (i=0; i<nphih; ++i)
12035 { dihedrals[i].atom1 = amber_data->
DihHAt1[i] / 3;
12036 dihedrals[i].atom2 = amber_data->
DihHAt2[i] / 3;
12037 dihedrals[i].atom3 = amber_data->
DihHAt3[i] / 3;
12038 dihedrals[i].atom4 = amber_data->
DihHAt4[i] / 3;
12039 dihedrals[i].dihedral_type = amber_data->
DihHNum[i] - 1;
12042 for (i=nphih; i<numDihedrals; ++i)
12043 { dihedrals[i].atom1 = amber_data->
DihAt1[i-nphih] / 3;
12044 dihedrals[i].atom2 = amber_data->
DihAt2[i-nphih] / 3;
12045 dihedrals[i].atom3 = amber_data->
DihAt3[i-nphih] / 3;
12046 dihedrals[i].atom4 = amber_data->
DihAt4[i-nphih] / 3;
12047 dihedrals[i].dihedral_type = amber_data->
DihNum[i-nphih] - 1;
12062 for (i=0; i<numDihedrals; ++i)
12063 {
if (dihedrals[i].atom3 < 0 || dihedrals[i].atom4 < 0)
12064 { dihedrals[i].atom3 = abs(dihedrals[i].atom3);
12065 dihedrals[i].atom4 = abs(dihedrals[i].atom4);
12068 {
if (dihedrals[i].atom1 < dihedrals[i].atom4)
12069 a1=dihedrals[i].atom1, a2=dihedrals[i].atom4;
12071 a1=dihedrals[i].atom4, a2=dihedrals[i].atom1;
12075 min=0, max=numExclusions-1;
12076 while (!found && min<=max)
12077 { index = (min+max)/2;
12078 if (exclusions[index].atom1 == a1)
12080 else if (exclusions[index].atom1 < a1)
12086 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
12089 for (j=index-1; j>=0 && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; --j);
12090 if (j<0 || exclusions[j].atom1!=a1)
12091 for (j=index; j<numExclusions && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; ++j);
12092 if (j<numExclusions && exclusions[j].atom1==a1)
12093 exclusions[j].modified = 1;
12095 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
12097 if (dihedrals[i].atom1>=numAtoms || dihedrals[i].atom2>=numAtoms ||
12098 dihedrals[i].atom3>=numAtoms || dihedrals[i].atom4>=numAtoms ||
12099 dihedrals[i].dihedral_type>=amber_data->
Nptra)
12100 {
char err_msg[128];
12101 sprintf(err_msg,
"DIHEDRAL # %d OVERFLOW IN PARM FILE", i+1);
12107 numRealBonds = numBonds;
12108 build_atom_status();
12127 read_parm(gromacsTopFile);
12129 #ifndef MEM_OPT_VERSION 12132 assignLCPOTypes( 3 );
12150 #ifdef MEM_OPT_VERSION 12151 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
12159 atoms =
new Atom[numAtoms];
12166 if (atoms == NULL || atomNames == NULL )
12167 NAMD_die(
"memory allocation failed when reading atom information");
12171 for (i=0; i<numAtoms; ++i) {
12172 char *resname = nameArena->getNewArray(11);
12173 char *atomname = nameArena->getNewArray(11);
12174 char *atomtype = nameArena->getNewArray(11);
12175 int resnum,typenum;
12178 if (resname == NULL || atomname == NULL || atomtype == NULL)
12179 NAMD_die(
"memory allocation failed when reading atom information");
12182 gf->
getAtom(i,&resnum,resname,atomname,atomtype,&typenum,
12185 atomNames[i].resname = resname;
12186 atomNames[i].atomname = atomname;
12187 atomNames[i].atomtype = atomtype;
12188 atoms[i].mass = mass;
12189 atoms[i].charge = charge;
12190 atoms[i].vdw_type = typenum;
12193 if ( tmpResLookup ) tmpResLookup =
12194 tmpResLookup->
append(
"MAIN", resnum+1, i);
12196 if(atomSegResids) {
12198 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
12199 one->
resid = resnum+1;
12208 }
else if (atoms[i].mass <= 0.05) {
12210 }
else if (atoms[i].mass < 1.0) {
12212 }
else if (atoms[i].mass <=3.5) {
12214 }
else if ((atomNames[i].atomname[0] ==
'O') &&
12215 (atoms[i].mass >= 14.0) &&
12216 (atoms[i].mass <= 18.0)) {
12223 bonds =
new Bond[numBonds];
12225 NAMD_die(
"memory allocation failed when reading bond information");
12226 for(i=0;i<numBonds;i++) {
12229 gf->
getBond(i,&atom1,&atom2,&type);
12230 bonds[i].atom1 = atom1;
12231 bonds[i].atom2 = atom2;
12232 bonds[i].bond_type = (
Index)type;
12237 angles =
new Angle[numAngles];
12238 if (angles == NULL)
12239 NAMD_die(
"memory allocation failed when reading angle information");
12240 for(i=0;i<numAngles;i++) {
12242 int atom1,atom2,atom3;
12243 gf->
getAngle(i,&atom1,&atom2,&atom3,&type);
12245 angles[i].atom1 = atom1;
12246 angles[i].atom2 = atom2;
12247 angles[i].atom3 = atom3;
12249 angles[i].angle_type=type;
12253 exclusions =
new Exclusion[numExclusions];
12317 dihedrals =
new Dihedral[numDihedrals];
12318 if (dihedrals == NULL)
12319 NAMD_die(
"memory allocation failed when reading dihedral information");
12320 for(i=0;i<numDihedrals;i++) {
12322 int atom1,atom2,atom3,atom4;
12323 gf->
getDihedral(i,&atom1,&atom2,&atom3,&atom4,&type);
12324 dihedrals[i].atom1 = atom1;
12325 dihedrals[i].atom2 = atom2;
12326 dihedrals[i].atom3 = atom3;
12327 dihedrals[i].atom4 = atom4;
12328 dihedrals[i].dihedral_type = type;
12341 const_cast<GromacsTopFile*
>(gf)->getPairLJArrays2(indxLJA, indxLJB, pairC6, pairC12);
12344 gromacsPair_type[i] = i;
12353 pointerToLJBeg =
new int[numAtoms];
12354 pointerToLJEnd =
new int[numAtoms];
12356 for(
int i=0; i < numAtoms; i++) {
12357 pointerToLJBeg[i] = -1;
12358 pointerToLJEnd[i] = -2;
12361 if(pointerToLJBeg[indxLJA[i]] == -1) {
12362 pointerToLJBeg[indxLJA[i]] = i;
12363 oldIndex = indxLJA[i];
12365 pointerToLJEnd[oldIndex] = i;
12378 const_cast<GromacsTopFile*
>(gf)->getPairGaussArrays2(indxGaussA, indxGaussB, gA, gMu1, giSigma1, gMu2, giSigma2, gRepulsive);
12381 pointerToGaussBeg =
new int[numAtoms];
12382 pointerToGaussEnd =
new int[numAtoms];
12383 for(
int i=0; i < numAtoms; i++) {
12384 pointerToGaussBeg[i] = -1;
12385 pointerToGaussEnd[i] = -2;
12389 if(pointerToGaussBeg[indxGaussA[i]] == -1) {
12390 pointerToGaussBeg[indxGaussA[i]] = i;
12391 oldIndex = indxGaussA[i];
12393 pointerToGaussEnd[oldIndex] = i;
12396 iout <<
iINFO <<
"Finished reading explicit pair from Gromacs file:\n" <<
12397 iINFO <<
"Found a total of: " << numPair <<
" explicit pairs--of which: " <<
12399 " are Gaussian style pairs.\n" <<
endi;
12403 #if GROMACS_EXCLUSIONS 12406 int* atom1 =
new int[numExclusions];
12407 int* atom2 =
new int[numExclusions];
12408 for(
int j=0; j<numExclusions;j++) {
12414 read_exclusions(atom1,atom2,numExclusions);
12478 numRealBonds = numBonds;
12479 build_atom_status();
12484 #ifndef MEM_OPT_VERSION 12498 #ifdef MEM_OPT_VERSION 12500 Index Molecule::insert_new_mass(
Real newMass){
12502 for(
int i=massPoolSize-1; i>=0; i--){
12503 if(fabs(atomMassPool[i]-newMass)<=1e-6)
12507 Real *tmp =
new Real[massPoolSize+1];
12508 tmp[massPoolSize] = newMass;
12509 memcpy((
void *)tmp, (
const void *)atomMassPool,
sizeof(
Real)*massPoolSize);
12510 delete [] atomMassPool;
12511 atomMassPool = tmp;
12513 return (
Index)(massPoolSize-1);
12516 void Molecule::addNewExclSigPool(
const vector<ExclusionSignature>& newExclSigPool){
12518 for(
int i=0; i<exclSigPoolSize; i++)
12519 tmpExclSigPool[i] = exclSigPool[i];
12520 for(
int i=0; i<newExclSigPool.size(); i++)
12521 tmpExclSigPool[i+exclSigPoolSize] = newExclSigPool[i];
12523 exclSigPoolSize += newExclSigPool.size();
12524 exclSigPool = tmpExclSigPool;
12528 msg->
put((
short)tupleType);
12529 msg->
put(numOffset);
12530 msg->
put(numOffset, offset);
12531 msg->
put(tupleParamType);
12540 msg->
get(numOffset);
12542 offset =
new int[numOffset];
12543 msg->
get(numOffset*
sizeof(
int), (
char *)offset);
12545 msg->
get(tupleParamType);
12551 for(
int i=0; i<bondCnt; i++)
12552 bondSigs[i].pack(msg);
12554 msg->
put(angleCnt);
12555 for(
int i=0; i<angleCnt; i++)
12556 angleSigs[i].pack(msg);
12558 msg->
put(dihedralCnt);
12559 for(
int i=0; i<dihedralCnt; i++)
12560 dihedralSigs[i].pack(msg);
12562 msg->
put(improperCnt);
12563 for(
int i=0; i<improperCnt; i++)
12564 improperSigs[i].pack(msg);
12566 msg->
put(crosstermCnt);
12567 for(
int i=0; i<crosstermCnt; i++)
12568 crosstermSigs[i].pack(msg);
12571 msg->
put(gromacsPairCnt);
12572 for(
int i=0; i<gromacsPairCnt; i++)
12573 gromacsPairSigs[i].pack(msg);
12578 delete [] bondSigs;
12581 for(
int i=0; i<bondCnt; i++)
12582 bondSigs[i].unpack(msg);
12583 }
else bondSigs = NULL;
12585 msg->
get(angleCnt);
12586 delete [] angleSigs;
12589 for(
int i=0; i<angleCnt; i++)
12590 angleSigs[i].unpack(msg);
12591 }
else angleSigs = NULL;
12593 msg->
get(dihedralCnt);
12594 delete [] dihedralSigs;
12597 for(
int i=0; i<dihedralCnt; i++)
12598 dihedralSigs[i].unpack(msg);
12599 }
else dihedralSigs = NULL;
12601 msg->
get(improperCnt);
12602 delete [] improperSigs;
12605 for(
int i=0; i<improperCnt; i++)
12606 improperSigs[i].unpack(msg);
12607 }
else improperSigs = NULL;
12609 msg->
get(crosstermCnt);
12610 delete [] crosstermSigs;
12611 if(crosstermCnt>0){
12613 for(
int i=0; i<crosstermCnt; i++)
12614 crosstermSigs[i].unpack(msg);
12615 }
else crosstermSigs = NULL;
12619 msg->
get(gromacsPairCnt);
12620 delete [] gromacsPairSigs;
12621 if(gromacsPairCnt>0){
12623 for(
int i=0; i<gromacsPairCnt; i++)
12624 gromacsPairSigs[i].unpack(msg);
12625 }
else gromacsPairSigs = NULL;
12639 origTupleCnt = bondCnt;
12640 tupleSigs= bondSigs;
12641 for(
int i=0; i<origTupleCnt; i++){
12646 delete [] tupleSigs;
12648 }
else if(bondCnt!=origTupleCnt){
12651 for(
int i=0; i<origTupleCnt; i++){
12653 newTupleSigs[idx] = tupleSigs[i];
12657 delete [] tupleSigs;
12658 bondSigs = newTupleSigs;
12664 origTupleCnt = angleCnt;
12665 tupleSigs = angleSigs;
12666 for(
int i=0; i<origTupleCnt; i++){
12671 delete [] tupleSigs;
12673 }
else if(angleCnt!=origTupleCnt){
12676 for(
int i=0; i<origTupleCnt; i++){
12678 newTupleSigs[idx] = tupleSigs[i];
12682 delete [] tupleSigs;
12683 angleSigs = newTupleSigs;
12689 origTupleCnt = dihedralCnt;
12690 tupleSigs = dihedralSigs;
12691 for(
int i=0; i<origTupleCnt; i++){
12695 if(dihedralCnt==0){
12696 delete [] tupleSigs;
12697 dihedralSigs = NULL;
12698 }
else if(dihedralCnt!=origTupleCnt){
12701 for(
int i=0; i<origTupleCnt; i++){
12703 newTupleSigs[idx] = tupleSigs[i];
12707 delete [] tupleSigs;
12708 dihedralSigs = newTupleSigs;
12715 origTupleCnt = improperCnt;
12716 tupleSigs = improperSigs;
12717 for(
int i=0; i<origTupleCnt; i++){
12721 if(improperCnt==0){
12722 delete [] tupleSigs;
12723 improperSigs = NULL;
12724 }
else if(improperCnt!=origTupleCnt){
12727 for(
int i=0; i<origTupleCnt; i++){
12729 newTupleSigs[idx] = tupleSigs[i];
12733 delete [] tupleSigs;
12734 improperSigs = newTupleSigs;
12740 origTupleCnt = crosstermCnt;
12741 tupleSigs = crosstermSigs;
12742 for(
int i=0; i<origTupleCnt; i++){
12746 if(crosstermCnt==0){
12747 delete [] tupleSigs;
12748 crosstermSigs = NULL;
12749 }
else if(crosstermCnt!=origTupleCnt){
12752 for(
int i=0; i<origTupleCnt; i++){
12754 newTupleSigs[idx] = tupleSigs[i];
12758 delete [] tupleSigs;
12759 crosstermSigs = newTupleSigs;
12766 origTupleCnt = gromacsPairCnt;
12767 tupleSigs = gromacsPairSigs;
12768 for(
int i=0; i<origTupleCnt; i++){
12772 if(gromacsPairCnt==0){
12773 delete [] tupleSigs;
12774 gromacsPairSigs = NULL;
12775 }
else if(gromacsPairCnt!=origTupleCnt){
12778 for(
int i=0; i<origTupleCnt; i++){
12780 newTupleSigs[idx] = tupleSigs[i];
12784 delete [] tupleSigs;
12785 gromacsPairSigs = newTupleSigs;
12795 for(
int i=0; i<fullExclCnt; i++){
12796 if(fullOffset[i]==0)
continue;
12801 delete [] fullOffset;
12803 }
else if(newCnt!=fullExclCnt){
12804 int *tmpOffset =
new int[newCnt];
12806 for(
int i=0; i<fullExclCnt; i++){
12807 if(fullOffset[i]==0)
continue;
12808 tmpOffset[newCnt] = fullOffset[i];
12811 delete [] fullOffset;
12812 fullOffset = tmpOffset;
12813 fullExclCnt = newCnt;
12818 for(
int i=0; i<modExclCnt; i++){
12819 if(modOffset[i]==0)
continue;
12824 delete [] modOffset;
12826 }
else if(newCnt!=modExclCnt){
12827 int *tmpOffset =
new int[newCnt];
12829 for(
int i=0; i<modExclCnt; i++){
12830 if(modOffset[i]==0)
continue;
12831 tmpOffset[newCnt] = modOffset[i];
12834 delete [] modOffset;
12835 modOffset = tmpOffset;
12836 modExclCnt = newCnt;
12851 int high = fullExclCnt-1;
12852 int mid = (low+high)/2;
12854 if(offset<fullOffset[mid]){
12856 mid = (high+low)/2;
12857 }
else if(offset>fullOffset[mid]){
12859 mid = (high+low)/2;
12865 if(retidx!=-1)
return retidx;
12869 high = modExclCnt-1;
12870 mid = (low+high)/2;
12872 if(offset<modOffset[mid]){
12874 mid = (high+low)/2;
12875 }
else if(offset>modOffset[mid]){
12877 mid = (high+low)/2;
12887 msg->
put(fullExclCnt);
12888 msg->
put(fullExclCnt, fullOffset);
12889 msg->
put(modExclCnt);
12890 msg->
put(modExclCnt, modOffset);
12894 msg->
get(fullExclCnt);
12895 delete [] fullOffset;
12896 fullOffset =
new int[fullExclCnt];
12897 msg->
get(fullExclCnt*
sizeof(
int), (
char *)fullOffset);
12898 msg->
get(modExclCnt);
12899 delete [] modOffset;
12900 modOffset =
new int[modExclCnt];
12901 msg->
get(modExclCnt*
sizeof(
int), (
char *)modOffset);
12902 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 12908 #endif // MOLECULE2_C defined = second object file
void build_gridforce_params(StringList *, StringList *, StringList *, StringList *, PDB *, char *)
void assign_improper_index(const char *, const char *, const char *, const char *, Improper *, int)
#define NAMD_EVENT_STOP(eon, id)
std::ostream & iINFO(std::ostream &s)
int get_mother_atom(int) const
void getAtom(int num, int *residue_number, char *residue_name, char *atom_name, char *atom_type, int *atom_typenum, Real *charge, Real *mass) const
void add_dcd_selection_file(int dcdIndex, char *userDcdFile)
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 *)
int checkexcl(int atom1, int atom2) const
void build_extra_bonds(Parameters *parameters, StringList *file)
void setOccupancyData(molfile_atom_t *atomarray)
int get_residue_size(const char *segid, int resid) const
void read_alch_unpert_dihedrals(FILE *)
vector< string > ResNames
Bool is_hydrogenGroupParent(int)
static GridforceGrid * new_grid(int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
BigReal temperaturefactor(void)
int lookup(const char *segid, int resid, int *begin, int *end) const
void send_Molecule(MOStream *)
void unpack(MIStream *msg)
static void pack_grid(GridforceGrid *grid, MOStream *msg)
void getBond(int num, int *atomi, int *atomj, int *bondtype) const
std::ostream & endi(std::ostream &s)
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
int add(const Elem &elem)
int getLCPOTypeAmber(char atomType[11], int numBonds)
void assign_bond_index(const char *, const char *, Bond *, bool *bond_found=nullptr)
void unpack(MIStream *msg)
#define COMPRESSED_PSF_VER
uint16_t find_or_create_dcd_selection_index(const char *keystr)
FourBodyConsts values[MAX_MULTIPLICITY]
void build_langevin_params(BigReal coupling, BigReal drudeCoupling, Bool doHydrogen)
static Units next(Units u)
void get_angle_params(Real *k, Real *theta0, Real *k_ub, Real *r_ub, Index index)
#define NAMD_EVENT_START(eon, id)
BigReal getEnergyTailCorr(const BigReal, const int)
int get_atom_from_name(const char *segid, int resid, const char *aname) const
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
void compute_LJcorrection_alternative()
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 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 getAngle(int num, int *atomi, int *atomj, int *atomk, int *angletype) const
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
NbtholePairValue * nbthole_array
int get_num_vdw_params(void)
FourBodyConsts values[MAX_MULTIPLICITY]
void build_fep_flags(StringList *, StringList *, PDB *, char *, const char *)
vector< string > AtomNames
ResidueLookupElem * append(const char *segid, int resid, int aid)
int sizeColumn(const vector< int > &v1, const vector< int > &v2)
void build_alch_unpert_bond_lists(char *)
int32 status
Atom status bit fields defined in structures.h.
int getNumDihedrals() const
struct OutputAtomRecord::floatVals fSet
int getNumGaussPair() const
void build_fixed_atoms(StringList *, StringList *, PDB *, char *)
void compute_LJcorrection()
int atomsInMigrationGroup
void print_atoms(Parameters *)
static GridforceGrid * unpack_grid(int gridnum, MIStream *msg)
struct OutputAtomRecord::integerVals iSet
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)
HashPool< HashString > atomTypePool
void build_exPressure_atoms(StringList *, StringList *, PDB *, char *)
void build_dcd_selection_list_pdb(int dcdIndex, char *userDcdInputFile)
void getDihedral(int num, int *atomi, int *atomj, int *atomk, int *atoml, int *type) const
void build_constant_forces(char *)
void unpack(MIStream *msg)
int get_atom_from_index_in_residue(const char *segid, int resid, int index) const
void assign_angle_index(const char *, const char *, const char *, Angle *, int)
MOStream * put(char data)
uint16_t find_dcd_selection_index(const char *keystr)
std::ostream & iERROR(std::ostream &s)
void get_bond_params(Real *k, Real *x0, Index index)
StringList * find(const char *name) const
void receive_Molecule(MIStream *)
void build_movdrag_params(StringList *, StringList *, StringList *, PDB *, char *)
struct OutputAtomRecord::shortVals sSet
void add_dcd_selection_freq(int dcdIndex, int freq)
void get_vdw_params(Real *sigma, Real *epsilon, Real *sigma14, Real *epsilon14, Index index)
int getNumExclusions() const
int getLCPOTypeCharmm(char atomType[11], int numBonds)
ResizeArray< int > atomIndex
Molecule(SimParameters *, Parameters *param)
void parse_dcd_selection_params(ConfigList *configList)
HashPool< AtomSigInfo > atomSigPool
void read_alch_unpert_angles(FILE *)