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;
4112 if (numOneFourNbTholes > 0 &&
simParams->alchOn) {
4113 NAMD_die(
"1-4 NbThole is incompatible with FEP/TI.\n");
4118 delete [] byAtomSize; byAtomSize = 0;
4119 delete [] byAtomSize2; byAtomSize2 = 0;
4125 for (i=0; i<numAtoms; i++)
4127 all_exclusions[i].
min = numAtoms;
4128 all_exclusions[i].max = -1;
4130 for (i=0; i<numTotalExclusions; i++)
4133 int a1 = exclusions[i].atom1;
4134 int a2 = exclusions[i].atom2;
4135 if ( numFixedAtoms && fixedAtomFlags[a1]
4136 && fixedAtomFlags[a2] )
continue;
4137 if ( all_exclusions[a1].min > a2 ) all_exclusions[a1].min = a2;
4138 if ( all_exclusions[a2].min > a1 ) all_exclusions[a2].min = a1;
4139 if ( a2 > all_exclusions[a1].max ) all_exclusions[a1].max = a2;
4140 if ( a1 > all_exclusions[a2].max ) all_exclusions[a2].max = a1;
4144 int maxDenseAllFull = 0;
4145 int numDenseAllFull = 0;
4146 for (i=0; i<numAtoms; i++) {
4147 int iInMiddle = ( i < all_exclusions[i].max &&
4148 i > all_exclusions[i].min ) ? 1 : 0;
4149 int s = all_exclusions[i].max - all_exclusions[i].min + 1;
4150 if ( s == fullExclusionsByAtom[i][0] + iInMiddle ) {
4151 if ( s > maxDenseAllFull ) maxDenseAllFull = s;
4152 all_exclusions[i].flags = (
char*)-1;
4154 all_exclusions[i].flags = 0;
4157 char *denseFullArray = exclArena->getNewArray(maxDenseAllFull);
4158 for ( i=0; i<maxDenseAllFull; ++i ) denseFullArray[i] =
EXCHCK_FULL;
4160 int exclmem = maxDenseAllFull;
4161 int maxExclusionFlags =
simParams->maxExclusionFlags;
4162 for (i=0; i<numAtoms; i++) {
4163 int s = all_exclusions[i].max - all_exclusions[i].min + 1;
4164 if ( all_exclusions[i].max != -1 ) {
4165 if ( all_exclusions[i].flags ) {
4166 all_exclusions[i].flags = denseFullArray;
4168 }
else if ( s < maxExclusionFlags ) {
4169 char *f = all_exclusions[i].flags = exclArena->getNewArray(s);
4170 for (
int k=0; k<s; ++k ) f[k] = 0;
4173 all_exclusions[i].flags = 0;
4176 all_exclusions[i].flags = (
char*)-1;
4180 iout <<
iINFO << numTotalExclusions <<
" exclusions consume " 4181 << exclmem <<
" bytes.\n" <<
endi;
4183 <<
" atoms sharing one array.\n" <<
endi;
4185 for (i=0; i<numTotalExclusions; i++)
4187 int a1 = exclusions[i].atom1;
4188 int a2 = exclusions[i].atom2;
4189 if ( numFixedAtoms && fixedAtomFlags[a1]
4190 && fixedAtomFlags[a2] )
continue;
4191 if ( exclusions[i].modified ) {
4192 if ( all_exclusions[a1].flags )
4193 all_exclusions[a1].flags[a2-all_exclusions[a1].min] =
EXCHCK_MOD;
4194 if ( all_exclusions[a2].flags )
4195 all_exclusions[a2].flags[a1-all_exclusions[a2].min] =
EXCHCK_MOD;
4197 if ( all_exclusions[a1].flags )
4198 all_exclusions[a1].flags[a2-all_exclusions[a1].min] =
EXCHCK_FULL;
4199 if ( all_exclusions[a2].flags )
4200 all_exclusions[a2].flags[a1-all_exclusions[a2].min] =
EXCHCK_FULL;
4229 void Molecule::build_exclusions()
4237 for (i=0; i<numExclusions; i++)
4239 exclusionSet.add(exclusions[i]);
4247 switch (exclude_flag)
4274 if (is_lonepairs_psf || is_drude_psf) {
4275 build_inherited_excl(
SCALED14 == exclude_flag);
4294 void Molecule::build_inherited_excl(
int modified) {
4296 int32 *bond1, *bond2, *bond3, *bond4, *bond5;
4297 int32 i, j, mid1, mid2, mid3, mid4;
4301 for (i = 0; i < numAtoms; i++) {
4303 if (!is_drude(i) && !is_lp(i))
continue;
4307 bond1 = bondsWithAtom[i];
4311 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4312 sprintf(err_msg,
"FOUND ISOLATED %s PARTICLE %d", idescrip, i+1);
4315 if (-1 != *(bond1+1)) {
4317 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4318 sprintf(err_msg,
"FOUND MULTIPLY LINKED %s PARTICLE %d",
4324 mid1 = bonds[*bond1].atom1;
4325 if (mid1 == i) mid1 = bonds[*bond1].atom2;
4328 if (is_drude(mid1) || is_lp(mid1) || is_hydrogen(mid1)) {
4330 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4331 sprintf(err_msg,
"PARENT ATOM %d of %s PARTICLE %d " 4332 "IS NOT HEAVY ATOM", mid1+1, idescrip, i+1);
4336 if (exclude_flag ==
NONE) {
4346 bond2 = bondsWithAtom[mid1];
4347 while (*bond2 != -1) {
4348 j = bonds[*bond2].atom1;
4349 if ((is_drude(j) || is_lp(j)) && j != mid1) {
4350 if (i < j) exclusionSet.add(
Exclusion(i, j));
4351 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4353 j = bonds[*bond2].atom2;
4354 if ((is_drude(j) || is_lp(j)) && j != mid1) {
4355 if (i < j) exclusionSet.add(
Exclusion(i, j));
4356 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4364 bond2 = bondsWithAtom[mid1];
4367 while (*bond2 != -1) {
4368 if (bonds[*bond2].atom1 == mid1) {
4369 mid2 = bonds[*bond2].atom2;
4372 mid2 = bonds[*bond2].atom1;
4382 if (exclude_flag ==
ONETWO) {
4392 bond3 = bondsWithAtom[mid2];
4393 while (*bond3 != -1) {
4394 j = bonds[*bond3].atom1;
4395 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4396 if (i < j) exclusionSet.add(
Exclusion(i, j));
4397 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4399 j = bonds[*bond3].atom2;
4400 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4401 if (i < j) exclusionSet.add(
Exclusion(i, j));
4402 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4410 bond3 = bondsWithAtom[mid2];
4413 while (*bond3 != -1) {
4415 if (bonds[*bond3].atom1 == mid2) {
4416 mid3 = bonds[*bond3].atom2;
4419 mid3 = bonds[*bond3].atom1;
4433 else if (mid3 < i) {
4439 bond4 = bondsWithAtom[mid3];
4440 while (*bond4 != -1) {
4441 j = bonds[*bond4].atom1;
4442 if ((is_drude(j) || is_lp(j)) && j != mid3) {
4443 if (i < j) exclusionSet.add(
Exclusion(i, j));
4444 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4446 j = bonds[*bond4].atom2;
4447 if ((is_drude(j) || is_lp(j)) && j != mid3) {
4448 if (i < j) exclusionSet.add(
Exclusion(i, j));
4449 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4457 bond4 = bondsWithAtom[mid3];
4460 while (*bond4 != -1) {
4462 if (bonds[*bond4].atom1 == mid3) {
4463 mid4 = bonds[*bond4].atom2;
4466 mid4 = bonds[*bond4].atom1;
4476 if (is_drude(mid4) || is_lp(mid4)) {
4481 else if (mid4 < i) {
4492 int modi = modified;
4494 int amin = (mid1 < mid4 ? mid1 : mid4);
4495 int amax = (mid1 >= mid4 ? mid1 : mid4);
4507 exclusionSet.add(
Exclusion(i, mid4, modi));
4509 else if (mid4 < i) {
4510 exclusionSet.add(
Exclusion(mid4, i, modi));
4515 bond5 = bondsWithAtom[mid4];
4516 while (*bond5 != -1) {
4517 j = bonds[*bond5].atom1;
4518 if ((is_drude(j) || is_lp(j)) && j != mid4) {
4519 if (i<j) exclusionSet.add(
Exclusion(i,j,modi));
4520 else if (j<i) exclusionSet.add(
Exclusion(j,i,modi));
4522 j = bonds[*bond5].atom2;
4523 if ((is_drude(j) || is_lp(j)) && j != mid4) {
4524 if (i<j) exclusionSet.add(
Exclusion(i,j,modi));
4525 else if (j<i) exclusionSet.add(
Exclusion(j,i,modi));
4555 void Molecule::build12excl(
void)
4561 for (i=0; i<numAtoms; i++)
4563 current_val = bondsWithAtom[i];
4566 while (*current_val != -1)
4568 if (bonds[*current_val].atom1 == i)
4570 if (i<bonds[*current_val].atom2)
4572 exclusionSet.add(
Exclusion(i,bonds[*current_val].atom2));
4577 if (i<bonds[*current_val].atom1)
4579 exclusionSet.add(
Exclusion(i,bonds[*current_val].atom1));
4595 void Molecule::build13excl(
void)
4597 int32 *bond1, *bond2;
4603 for (i=0; i<numAtoms; i++)
4605 bond1 = bondsWithAtom[i];
4608 while (*bond1 != -1)
4610 if (bonds[*bond1].atom1 == i)
4612 middle_atom=bonds[*bond1].atom2;
4616 middle_atom=bonds[*bond1].atom1;
4619 bond2 = bondsWithAtom[middle_atom];
4623 while (*bond2 != -1)
4625 if (bonds[*bond2].atom1 == middle_atom)
4627 if (i < bonds[*bond2].atom2)
4629 exclusionSet.add(
Exclusion(i,bonds[*bond2].atom2));
4634 if (i < bonds[*bond2].atom1)
4636 exclusionSet.add(
Exclusion(i,bonds[*bond2].atom1));
4656 void Molecule::build14excl(
int modified)
4658 int32 *bond1, *bond2, *bond3;
4663 for (i=0; i<numAtoms; i++)
4665 if (is_drude(i) || is_lp(i))
continue;
4668 bond1 = bondsWithAtom[i];
4670 while (*bond1 != -1)
4672 if (bonds[*bond1].atom1 == i)
4674 mid1=bonds[*bond1].atom2;
4678 mid1=bonds[*bond1].atom1;
4681 bond2 = bondsWithAtom[mid1];
4684 while (*bond2 != -1)
4686 if (bonds[*bond2].atom1 == mid1)
4688 mid2 = bonds[*bond2].atom2;
4692 mid2 = bonds[*bond2].atom1;
4704 bond3=bondsWithAtom[mid2];
4707 while (*bond3 != -1)
4709 if (bonds[*bond3].atom1 == mid2)
4711 int j = bonds[*bond3].atom2;
4717 if (i < j && !is_drude(j) && !is_lp(j))
4719 exclusionSet.add(
Exclusion(i,j,modified));
4724 int j = bonds[*bond3].atom1;
4730 if (i < j && !is_drude(j) && !is_lp(j))
4732 exclusionSet.add(
Exclusion(i,j,modified));
4754 void Molecule::stripFepExcl(
void)
4760 for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
4762 int t1 = get_fep_type(exclIter->atom1);
4763 int t2 = get_fep_type(exclIter->atom2);
4764 if (t1 && t2 && t1 !=t2 && abs(t1-t2) != 2) {
4765 fepExclusionSet.
add(*exclIter);
4768 }
else if (
simParams->pairInteractionOn ) {
4769 for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
4771 int ifep_type = get_fep_type(exclIter->atom1);
4772 int jfep_type = get_fep_type(exclIter->atom2);
4775 if (ifep_type != 1 || jfep_type != 1) {
4776 fepExclusionSet.
add(*exclIter);
4781 if (!(ifep_type == 1 && jfep_type == 2) &&
4782 !(ifep_type == 2 && jfep_type == 1)) {
4783 fepExclusionSet.
add(*exclIter);
4790 for ( fepIter=fepIter.begin(); fepIter != fepIter.end(); fepIter++ )
4792 exclusionSet.del(*fepIter);
4806 if ( (psf_file =
Fopen(fname,
"r")) == NULL)
4809 sprintf(err_msg,
"UNABLE TO OPEN THE COMPRESSED .psf FILE %s", fname);
4818 NAMD_die(
"The compressed psf file format is incorrect, please re-generate!\n");
4820 float psfVer = 0.0f;
4821 sscanf(buffer,
"FORMAT VERSION: %f\n", &psfVer);
4823 NAMD_die(
"The compressed psf file format is incorrect, please re-generate!\n");
4828 NAMD_die(
"UNABLE TO FIND NSEGMENTNAMES");
4829 sscanf(buffer,
"%d", &segNamePoolSize);
4831 if(segNamePoolSize!=0)
4833 for(
int i=0; i<segNamePoolSize; i++){
4835 sscanf(buffer,
"%s", strBuf);
4836 segNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4840 for(
int i=0; i<segNamePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4845 NAMD_die(
"UNABLE TO FIND NRESIDUENAMES");
4846 sscanf(buffer,
"%d", &resNamePoolSize);
4848 if(resNamePoolSize!=0)
4850 for(
int i=0; i<resNamePoolSize; i++){
4852 sscanf(buffer,
"%s", strBuf);
4853 resNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4857 for(
int i=0; i<resNamePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4862 NAMD_die(
"UNABLE TO FIND NATOMNAMES");
4863 sscanf(buffer,
"%d", &atomNamePoolSize);
4864 if(atomNamePoolSize!=0)
4866 for(
int i=0; i<atomNamePoolSize; i++){
4868 sscanf(buffer,
"%s", strBuf);
4869 atomNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4875 NAMD_die(
"UNABLE TO FIND NATOMTYPES");
4876 sscanf(buffer,
"%d", &atomTypePoolSize);
4878 if(atomTypePoolSize!=0)
4880 for(
int i=0; i<atomTypePoolSize; i++){
4882 sscanf(buffer,
"%s", strBuf);
4883 atomTypePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4887 for(
int i=0; i<atomTypePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4892 NAMD_die(
"UNABLE TO FIND NCHARGES");
4893 sscanf(buffer,
"%d", &chargePoolSize);
4894 if(chargePoolSize!=0)
4895 atomChargePool =
new Real[chargePoolSize];
4896 for(
int i=0; i<chargePoolSize; i++){
4898 sscanf(buffer,
"%f", atomChargePool+i);
4903 NAMD_die(
"UNABLE TO FIND NMASSES");
4904 sscanf(buffer,
"%d", &massPoolSize);
4906 atomMassPool =
new Real[massPoolSize];
4907 for(
int i=0; i<massPoolSize; i++){
4909 sscanf(buffer,
"%f", atomMassPool+i);
4914 NAMD_die(
"UNABLE TO FIND ATOMSIGS");
4915 sscanf(buffer,
"%d", &atomSigPoolSize);
4918 int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
4921 for(
int i=0; i<atomSigPoolSize; i++){
4925 NAMD_die(
"UNABLE TO FIND NBONDSIGS");
4926 sscanf(buffer,
"%d", &typeCnt);
4931 for(
int j=0; j<typeCnt; j++){
4933 sscanf(buffer,
"%d | %d | %d", &tmp1, &ttype, &tisReal);
4935 oneSig.offset[0] = tmp1;
4937 if(tisReal) numRealBonds++;
4943 NAMD_die(
"UNABLE TO FIND NTHETASIGS");
4944 sscanf(buffer,
"%d", &typeCnt);
4949 for(
int j=0; j<typeCnt; j++){
4951 sscanf(buffer,
"%d %d | %d | %d", &tmp1, &tmp2, &ttype, &tisReal);
4953 oneSig.offset[0] = tmp1;
4954 oneSig.offset[1] = tmp2;
4960 NAMD_die(
"UNABLE TO FIND NPHISIGS");
4961 sscanf(buffer,
"%d", &typeCnt);
4966 for(
int j=0; j<typeCnt; j++){
4968 sscanf(buffer,
"%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
4970 oneSig.offset[0] = tmp1;
4971 oneSig.offset[1] = tmp2;
4972 oneSig.offset[2] = tmp3;
4978 NAMD_die(
"UNABLE TO FIND NIMPHISIGS");
4979 sscanf(buffer,
"%d", &typeCnt);
4984 for(
int j=0; j<typeCnt; j++){
4986 sscanf(buffer,
"%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
4988 oneSig.offset[0] = tmp1;
4989 oneSig.offset[1] = tmp2;
4990 oneSig.offset[2] = tmp3;
4996 NAMD_die(
"UNABLE TO FIND NCRTERMSIGS");
4997 sscanf(buffer,
"%d", &typeCnt);
5002 for(
int j=0; j<typeCnt; j++){
5004 sscanf(buffer,
"%d %d %d %d %d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &ttype, &tisReal);
5006 oneSig.offset[0] = tmp1;
5007 oneSig.offset[1] = tmp2;
5008 oneSig.offset[2] = tmp3;
5009 oneSig.offset[3] = tmp4;
5010 oneSig.offset[4] = tmp5;
5011 oneSig.offset[5] = tmp6;
5012 oneSig.offset[6] = tmp7;
5020 NAMD_die(
"UNABLE TO FIND NEXCLSIGS");
5022 sscanf(buffer,
"%d", &exclSigPoolSize);
5024 vector<int> fullExcls;
5025 vector<int> modExcls;
5026 for(
int i=0; i<exclSigPoolSize; i++){
5028 for(
int j=0; j<fullExclCnt; j++)
5031 for(
int j=0; j<modExclCnt; j++)
5035 exclSigPool[i].setOffsets(fullExcls, modExcls);
5044 NAMD_die(
"UNABLE TO FIND NCLUSTERS");
5046 sscanf(buffer,
"%d", &numClusters);
5051 sscanf(buffer,
"%d", &numAtoms);
5055 NAMD_die(
"UNABLE TO FIND NHYDROGENGROUP");
5056 sscanf(buffer,
"%d", &numHydrogenGroups);
5060 NAMD_die(
"UNABLE TO FIND MAXHYDROGENGROUPSIZE");
5061 sscanf(buffer,
"%d", &maxHydrogenGroupSize);
5064 NAMD_die(
"UNABLE TO FIND NMIGRATIONGROUP");
5065 sscanf(buffer,
"%d", &numMigrationGroups);
5068 NAMD_die(
"UNABLE TO FIND MAXMIGRATIONGROUPSIZE");
5069 sscanf(buffer,
"%d", &maxMigrationGroupSize);
5071 int inputRigidType = -1;
5074 NAMD_die(
"UNABLE TO FIND RIGIDBONDTYPE");
5075 sscanf(buffer,
"%d", &inputRigidType);
5078 if(
simParams->rigidBonds != inputRigidType){
5079 const char *tmpstr[]={
"RIGID_NONE",
"RIGID_ALL",
"RIGID_WATER"};
5081 sprintf(errmsg,
"RIGIDBOND TYPE MISMATCH BETWEEN INPUT (%s) AND CURRENT RUN (%s)",
5082 tmpstr[inputRigidType], tmpstr[
simParams->rigidBonds]);
5090 NAMD_die(
"UNABLE TO FIND OCCUPANCYVALID");
5091 sscanf(buffer,
"%d", &isOccupancyValid);
5094 NAMD_die(
"UNABLE TO FIND TEMPFACTORVALID");
5095 sscanf(buffer,
"%d", &isBFactorValid);
5103 build_extra_bonds(params, cfgList->
find(
"extraBondsFile"));
5107 NAMD_die(
"UNABLE TO FIND DIHEDRALPARAMARRAY");
5116 NAMD_die(
"UNABLE TO FIND IMPROPERPARAMARRAY");
5133 void Molecule::read_binary_atom_info(
int fromAtomID,
int toAtomID,
InputAtomList& inAtoms){
5134 int numAtomsPar = toAtomID-fromAtomID+1;
5135 CmiAssert(numAtomsPar > 0);
5136 CmiAssert(inAtoms.
size() == numAtomsPar);
5155 eachAtomMass = NULL;
5156 eachAtomCharge = NULL;
5158 eachAtomExclSig = NULL;
5181 FILE *perAtomFile = fopen(
simParams->binAtomFile,
"rb");
5182 if (perAtomFile==NULL) {
5184 sprintf(err_msg,
"UNABLE TO OPEN THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s",
simParams->binAtomFile);
5190 flipNum((
char *)&rMagicNum,
sizeof(
int), 1);
5192 fread(&fMagicNum,
sizeof(
int), 1, perAtomFile);
5193 if (fMagicNum==magicNum) {
5195 }
else if (fMagicNum==rMagicNum) {
5199 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS CORRUPTED",
simParams->binAtomFile);
5203 float verNum = 0.0f;
5204 fread(&verNum,
sizeof(
float), 1, perAtomFile);
5205 if (needFlip)
flipNum((
char *)&verNum,
sizeof(
float), 1);
5208 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n",
simParams->binAtomFile);
5213 fread(&recSize,
sizeof(
int), 1, perAtomFile);
5214 if(needFlip)
flipNum((
char *)&recSize,
sizeof(
int), 1);
5217 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM RECORD SIZE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n",
simParams->binAtomFile);
5221 const int BUFELEMS = 32*1024;
5226 if ( _fseeki64(perAtomFile,startbyte,SEEK_CUR) )
5228 if ( fseeko(perAtomFile,startbyte,SEEK_CUR) )
5232 sprintf(errmsg,
"Error on seeking binary file %s",
simParams->binAtomFile);
5238 int atomsCnt = numAtomsPar;
5241 while(atomsCnt >= BUFELEMS) {
5242 if ( fread((
char *)elemsBuf,
sizeof(
OutputAtomRecord), BUFELEMS, perAtomFile) != BUFELEMS ) {
5244 sprintf(errmsg,
"Error on reading binary file %s",
simParams->binAtomFile);
5248 for(
int i=0; i<BUFELEMS; i++, curIdx++, oneRec++) {
5250 int aid = curIdx+fromAtomID;
5251 if(needFlip) oneRec->
flip();
5252 load_one_inputatom(aid, oneRec, fAtom);
5254 atomsCnt -= BUFELEMS;
5257 if ( fread(elemsBuf,
sizeof(
OutputAtomRecord), atomsCnt, perAtomFile) != atomsCnt ) {
5259 sprintf(errmsg,
"Error on reading binary file %s",
simParams->binAtomFile);
5263 for(
int i=curIdx; i<numAtomsPar; i++, oneRec++) {
5265 int aid = i+fromAtomID;
5266 if(needFlip) oneRec->
flip();
5267 load_one_inputatom(aid,oneRec,fAtom);
5270 if ( fclose(perAtomFile) ) {
5272 sprintf(errmsg,
"Error on closing binary file %s",
simParams->binAtomFile);
5281 is_atom_fixed(fromAtomID, &listIdx);
5282 for(
int i=listIdx; i<fixedAtomsSet->size(); i++){
5283 const AtomSet one = fixedAtomsSet->item(i);
5285 int sAtomId = one.aid1>fromAtomID ? one.aid1:fromAtomID;
5286 int eAtomId = one.aid2>toAtomID? toAtomID:one.aid2;
5287 for(
int j=sAtomId; j<=eAtomId; j++)
5288 inAtoms[j-fromAtomID].atomFixed = 1;
5295 char *thisAtomName = NULL;
5357 }
else if (thisAtomMass <= 0.05) {
5359 }
else if (thisAtomMass < 1.0) {
5361 }
else if (thisAtomMass <= 3.5) {
5363 }
else if (thisAtomName[0]==
'O' &&
5364 (thisAtomMass >= 14.0) && (thisAtomMass <= 18.0)) {
5378 void Molecule::build_excl_check_signatures(){
5380 for(
int i=0; i<exclSigPoolSize; i++){
5388 int fullMin, fullMax, modMin, modMax;
5396 if(fullMin < modMin)
5397 sigChk->
min = fullMin;
5399 sigChk->
min = modMin;
5400 if(fullMax < modMax)
5401 sigChk->
max = modMax;
5403 sigChk->
max = fullMax;
5411 iout <<
iWARN <<
"an empty exclusion signature with index " 5412 << i <<
"!\n" <<
endi;
5417 sigChk->
flags =
new char[sigChk->
max-sigChk->
min+1];
5418 memset(sigChk->
flags, 0,
sizeof(
char)*(sigChk->
max-sigChk->
min+1));
5439 void Molecule::load_atom_set(
StringList *setfile,
const char *setname,
5440 int *numAtomsInSet, AtomSetList **atomsSet)
const {
5441 if(setfile == NULL) {
5443 sprintf(errmsg,
"The text input file for %s atoms is not found!", setname);
5446 FILE *ifp = fopen(setfile->
data,
"r");
5450 sprintf(errmsg,
"ERROR IN OPENING %s ATOMS FILE: %s\n", setname, setfile->
data);
5455 int numLocalAtoms = 0;
5457 AtomSetList *localAtomsSet =
new AtomSetList();
5462 bool hasDash =
false;
5463 for(
int i=0; oneline[i] && i<128; i++){
5464 if(oneline[i]==
'-') {
5470 sscanf(oneline,
"%d-%d", &(one.aid1), &(one.aid2));
5471 if(one.aid1>one.aid2 || one.aid1<0 || one.aid2<0) {
5473 sprintf(errmsg,
"The input for %s atoms is wrong: %s\n", setname, oneline);
5476 numLocalAtoms += (one.aid2-one.aid1+1);
5478 sscanf(oneline,
"%d", &(one.aid1));
5481 sprintf(errmsg,
"The input for %s atoms is wrong: %s\n", setname, oneline);
5484 one.aid2 = one.aid1;
5487 localAtomsSet->add(one);
5491 std::sort(localAtomsSet->begin(), localAtomsSet->end());
5493 *numAtomsInSet = numLocalAtoms;
5494 *atomsSet = localAtomsSet;
5497 void Molecule::load_fixed_atoms(
StringList *fixedfile){
5498 load_atom_set(fixedfile,
"FIXED", &numFixedAtoms, &fixedAtomsSet);
5501 void Molecule::load_constrained_atoms(
StringList *constrainedfile){
5502 load_atom_set(constrainedfile,
"CONSTRAINED", &numConstraints, &constrainedAtomsSet);
5505 Bool Molecule::is_atom_in_set(AtomSetList *localAtomsSet,
int aid,
int *listIdx)
const {
5506 int idx = localAtomsSet->size();
5508 int lIdx = localAtomsSet->size()-1;
5510 while(rIdx <= lIdx){
5511 int mIdx = (rIdx+lIdx)/2;
5512 const AtomSet one = localAtomsSet->item(mIdx);
5518 }
else if(aid > one.aid1){
5522 if(listIdx) *listIdx = mIdx;
5529 if(listIdx) *listIdx = mIdx;
5535 if(listIdx) *listIdx = idx;
5553 #ifdef MEM_OPT_VERSION 5554 DebugM(3,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5563 <<
"******************************************\n" \
5564 <<
"NUM NAME TYPE RES MASS CHARGE CHARGE FEP-CHARGE" \
5565 <<
"SIGMA EPSILON SIGMA14 EPSILON14\n" \
5569 for (i=0; i<numAtoms; i++)
5571 const int vdw_type =
simParams->soluteScalingOn ?
5572 ((atoms[i].vdw_type >= LJtypecount) ?
5573 ss_vdw_type[atoms[i].vdw_type-LJtypecount] : atoms[i].vdw_type) : atoms[i].vdw_type;
5574 params->
get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, vdw_type);
5576 DebugM(3,i+1 <<
" " << atomNames[i].atomname \
5577 <<
" " << atomNames[i].atomtype <<
" " \
5578 << atomNames[i].resname <<
" " << atoms[i].mass \
5579 <<
" " << atoms[i].charge <<
" " << sigma \
5580 <<
" " << epsilon <<
" " << sigma14 \
5581 <<
" " << epsilon14 <<
"\n" \
5599 #ifdef MEM_OPT_VERSION 5600 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5606 DebugM(2,
"BOND LIST\n" <<
"********************************\n" \
5607 <<
"ATOM1 ATOM2 TYPE1 TYPE2 k x0" \
5610 for (i=0; i<numBonds; i++)
5614 DebugM(2,bonds[i].atom1+1 <<
" " \
5615 << bonds[i].atom2+1 <<
" " \
5616 << atomNames[bonds[i].atom1].atomtype <<
" " \
5617 << atomNames[bonds[i].atom2].atomtype <<
" " << k \
5618 <<
" " << x0 <<
endi);
5636 #ifdef MEM_OPT_VERSION 5637 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5641 DebugM(2,
"EXPLICIT EXCLUSION LIST\n" \
5642 <<
"********************************\n" \
5646 for (i=0; i<numExclusions; i++)
5648 DebugM(2,exclusions[i].atom1+1 <<
" " \
5649 << exclusions[i].atom2+1 <<
endi);
5666 #ifdef MEM_OPT_VERSION 5672 msg->
put(massPoolSize);
5673 msg->
put(massPoolSize, atomMassPool);
5675 msg->
put(chargePoolSize);
5676 msg->
put(chargePoolSize, atomChargePool);
5679 msg->
put(atomSigPoolSize);
5680 for(
int i=0; i<atomSigPoolSize; i++)
5684 msg->
put(exclSigPoolSize);
5685 for(
int i=0; i<exclSigPoolSize; i++)
5686 exclSigPool[i].pack(msg);
5688 msg->
put(numHydrogenGroups);
5689 msg->
put(maxHydrogenGroupSize);
5690 msg->
put(numMigrationGroups);
5691 msg->
put(maxMigrationGroupSize);
5692 msg->
put(isOccupancyValid);
5693 msg->
put(isBFactorValid);
5696 msg->
put(atomNamePoolSize);
5697 for(
int i=0; i<atomNamePoolSize;i++) {
5704 int numFixedAtomsSet = fixedAtomsSet->size();
5705 msg->
put(numFixedAtoms);
5706 msg->
put(numFixedAtomsSet);
5707 msg->
put(numFixedAtomsSet*
sizeof(AtomSet), (
char *)(fixedAtomsSet->begin()));
5711 int numConstrainedAtomsSet = constrainedAtomsSet->size();
5712 msg->
put(numConstraints);
5713 msg->
put(numConstrainedAtomsSet);
5714 msg->
put(numConstrainedAtomsSet*
sizeof(AtomSet), (
char *)(constrainedAtomsSet->begin()));
5719 msg->
put(numAtoms*
sizeof(
Atom), (
char*)atoms);
5722 msg->
put(numRealBonds);
5727 msg->
put(numBonds*
sizeof(
Bond), (
char*)bonds);
5731 msg->
put(numAngles);
5734 msg->
put(numAngles*
sizeof(
Angle), (
char*)angles);
5738 msg->
put(numDihedrals);
5741 msg->
put(numDihedrals*
sizeof(
Dihedral), (
char*)dihedrals);
5745 msg->
put(num_alch_unpert_Bonds);
5746 msg->
put(num_alch_unpert_Bonds*
sizeof(
Bond), (
char*)alch_unpert_bonds);
5748 msg->
put(num_alch_unpert_Angles);
5749 msg->
put(num_alch_unpert_Angles*
sizeof(
Angle), (
char*)alch_unpert_angles);
5751 msg->
put(num_alch_unpert_Dihedrals);
5752 msg->
put(num_alch_unpert_Dihedrals*
sizeof(
Dihedral), (
char*)alch_unpert_dihedrals);
5756 msg->
put(numImpropers);
5759 msg->
put(numImpropers*
sizeof(
Improper), (
char*)impropers);
5763 msg->
put(numCrossterms);
5766 msg->
put(numCrossterms*
sizeof(
Crossterm), (
char*)crossterms);
5770 msg->
put(numDonors);
5773 msg->
put(numDonors*
sizeof(
Bond), (
char*)donors);
5777 msg->
put(numAcceptors);
5780 msg->
put(numAcceptors*
sizeof(
Bond), (
char*)acceptors);
5784 msg->
put(numExclusions);
5787 msg->
put(numExclusions*
sizeof(
Exclusion), (
char*)exclusions);
5792 msg->
put(numConstraints);
5794 msg->
put(numAtoms, consIndexes);
5798 msg->
put(numConstraints*
sizeof(ConstraintParams), (
char*)consParams);
5807 DebugM(3,
"Sending gridforce info\n" <<
endi);
5808 msg->
put(numGridforceGrids);
5810 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
5811 msg->
put(numGridforces[gridnum]);
5812 msg->
put(numAtoms, gridfrcIndexes[gridnum]);
5813 if (numGridforces[gridnum])
5815 msg->
put(numGridforces[gridnum]*
sizeof(GridforceParams), (
char*)gridfrcParams[gridnum]);
5826 msg->
put(numStirredAtoms);
5828 msg->
put(numAtoms, stirIndexes);
5830 if (numStirredAtoms)
5833 msg->
put(numStirredAtoms*
sizeof(StirParams), (
char*)stirParams);
5840 msg->
put(numMovDrag);
5841 msg->
put(numAtoms, movDragIndexes);
5844 msg->
put(numMovDrag*
sizeof(MovDragParams), (
char*)movDragParams);
5850 msg->
put(numRotDrag);
5851 msg->
put(numAtoms, rotDragIndexes);
5854 msg->
put(numRotDrag*
sizeof(RotDragParams), (
char*)rotDragParams);
5860 msg->
put(numConsTorque);
5861 msg->
put(numAtoms, consTorqueIndexes);
5864 msg->
put(numConsTorque*
sizeof(ConsTorqueParams), (
char*)consTorqueParams);
5870 { msg->
put(numConsForce);
5871 msg->
put(numAtoms, consForceIndexes);
5873 msg->
put(numConsForce*
sizeof(
Vector), (
char*)consForce);
5878 msg->
put(numMolecules);
5879 msg->
put(numLargeMolecules);
5880 msg->
put(numAtoms, moleculeAtom);
5881 msg->
put(numMolecules+1, moleculeStartIndex);
5885 msg->
put(numExPressureAtoms);
5886 msg->
put(numAtoms, exPressureAtomFlags);
5889 #ifndef MEM_OPT_VERSION 5893 msg->
put(numAtoms, langevinParams);
5899 msg->
put(numFixedAtoms);
5900 msg->
put(numAtoms, fixedAtomFlags);
5901 msg->
put(numFixedRigidBonds);
5906 msg->
put(numAtoms, qmAtomGroup);
5907 msg->
put(qmNumQMAtoms);
5908 msg->
put(qmNumQMAtoms, qmAtmChrg);
5909 msg->
put(qmNumQMAtoms, qmAtmIndx);
5911 msg->
put(qmNumBonds);
5912 msg->
put(qmMeNumBonds);
5913 msg->
put(qmMeNumBonds, qmMeMMindx);
5914 msg->
put(qmMeNumBonds, qmMeQMGrp);
5916 msg->
put(qmNumGrps);
5917 msg->
put(qmNumGrps, qmGrpID);
5918 msg->
put(qmNumGrps, qmCustPCSizes);
5919 msg->
put(qmTotCustPCs);
5920 msg->
put(qmTotCustPCs, qmCustomPCIdxs);
5926 msg->
put(numFepInitial);
5927 msg->
put(numFepFinal);
5928 msg->
put(numAtoms*
sizeof(
char), (
char*)fepAtomFlags);
5933 msg->
put(numAtoms*
sizeof(
char), (
char*)ssAtomFlags);
5934 msg->
put(ss_num_vdw_params);
5936 msg->
put(numAtoms*
sizeof(
int), (
char*)ss_index);
5939 #ifdef OPENATOM_VERSION 5942 msg->
put(numFepInitial);
5943 msg->
put(numAtoms*
sizeof(
char), (
char*)fepAtomFlags);
5945 #endif //OPENATOM_VERSION 5948 msg->
put(is_lonepairs_psf);
5949 if (is_lonepairs_psf) {
5950 msg->
put(numLphosts);
5951 msg->
put(numLphosts*
sizeof(
Lphost), (
char*)lphosts);
5953 msg->
put(is_drude_psf);
5956 msg->
put(numTholes);
5957 msg->
put(numTholes*
sizeof(
Thole), (
char*)tholes);
5958 msg->
put(numAnisos);
5959 msg->
put(numAnisos*
sizeof(
Aniso), (
char*)anisos);
5961 msg->
put(numZeroMassAtoms);
5966 msg->
put(numAtoms, (
int*)lcpoParamType);
5977 msg->
put((numAtoms),pointerToLJBeg);
5978 msg->
put((numAtoms),pointerToLJEnd);
5988 msg->
put((numAtoms),pointerToGaussBeg);
5989 msg->
put((numAtoms),pointerToGaussEnd);
5997 #ifdef MEM_OPT_VERSION 5999 build_excl_check_signatures();
6002 numBonds = numCalcBonds = 0;
6003 numAngles = numCalcAngles = 0;
6004 numDihedrals = numCalcDihedrals = 0;
6005 numImpropers = numCalcImpropers = 0;
6006 numCrossterms = numCalcCrossterms = 0;
6007 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;
6015 build_lists_by_atom();
6035 #ifdef MEM_OPT_VERSION 6039 msg->
get(massPoolSize);
6040 if(atomMassPool)
delete [] atomMassPool;
6041 atomMassPool =
new Real[massPoolSize];
6042 msg->
get(massPoolSize, atomMassPool);
6044 msg->
get(chargePoolSize);
6045 if(atomChargePool)
delete [] atomChargePool;
6046 atomChargePool =
new Real[chargePoolSize];
6047 msg->
get(chargePoolSize, atomChargePool);
6050 msg->
get(atomSigPoolSize);
6053 for(
int i=0; i<atomSigPoolSize; i++)
6057 msg->
get(exclSigPoolSize);
6058 if(exclSigPool)
delete [] exclSigPool;
6060 for(
int i=0; i<exclSigPoolSize; i++)
6061 exclSigPool[i].unpack(msg);
6063 msg->
get(numHydrogenGroups);
6064 msg->
get(maxHydrogenGroupSize);
6065 msg->
get(numMigrationGroups);
6066 msg->
get(maxMigrationGroupSize);
6067 msg->
get(isOccupancyValid);
6068 msg->
get(isBFactorValid);
6071 msg->
get(atomNamePoolSize);
6073 for(
int i=0; i<atomNamePoolSize;i++) {
6081 int numFixedAtomsSet;
6082 msg->
get(numFixedAtoms);
6083 msg->
get(numFixedAtomsSet);
6084 fixedAtomsSet =
new AtomSetList(numFixedAtomsSet);
6085 msg->
get(numFixedAtomsSet*
sizeof(AtomSet), (
char *)(fixedAtomsSet->begin()));
6089 int numConstrainedAtomsSet;
6090 msg->
get(numConstraints);
6091 msg->
get(numConstrainedAtomsSet);
6092 constrainedAtomsSet =
new AtomSetList(numConstrainedAtomsSet);
6093 msg->
get(numConstrainedAtomsSet*
sizeof(AtomSet), (
char *)(constrainedAtomsSet->begin()));
6098 atoms=
new Atom[numAtoms];
6099 msg->
get(numAtoms*
sizeof(
Atom), (
char*)atoms);
6102 msg->
get(numRealBonds);
6107 bonds=
new Bond[numBonds];
6108 msg->
get(numBonds*
sizeof(
Bond), (
char*)bonds);
6112 msg->
get(numAngles);
6116 angles=
new Angle[numAngles];
6117 msg->
get(numAngles*
sizeof(
Angle), (
char*)angles);
6121 msg->
get(numDihedrals);
6124 delete [] dihedrals;
6125 dihedrals=
new Dihedral[numDihedrals];
6126 msg->
get(numDihedrals*
sizeof(
Dihedral), (
char*)dihedrals);
6130 msg->
get(num_alch_unpert_Bonds);
6131 alch_unpert_bonds=
new Bond[num_alch_unpert_Bonds];
6132 msg->
get(num_alch_unpert_Bonds*
sizeof(
Bond), (
char*)alch_unpert_bonds);
6134 msg->
get(num_alch_unpert_Angles);
6135 alch_unpert_angles=
new Angle[num_alch_unpert_Angles];
6136 msg->
get(num_alch_unpert_Angles*
sizeof(
Angle), (
char*)alch_unpert_angles);
6138 msg->
get(num_alch_unpert_Dihedrals);
6139 alch_unpert_dihedrals=
new Dihedral[num_alch_unpert_Dihedrals];
6140 msg->
get(num_alch_unpert_Dihedrals*
sizeof(
Dihedral), (
char*)alch_unpert_dihedrals);
6144 msg->
get(numImpropers);
6147 delete [] impropers;
6148 impropers=
new Improper[numImpropers];
6149 msg->
get(numImpropers*
sizeof(
Improper), (
char*)impropers);
6153 msg->
get(numCrossterms);
6156 delete [] crossterms;
6157 crossterms=
new Crossterm[numCrossterms];
6158 msg->
get(numCrossterms*
sizeof(
Crossterm), (
char*)crossterms);
6162 msg->
get(numDonors);
6166 donors=
new Bond[numDonors];
6167 msg->
get(numDonors*
sizeof(
Bond), (
char*)donors);
6171 msg->
get(numAcceptors);
6174 delete [] acceptors;
6175 acceptors=
new Bond[numAcceptors];
6176 msg->
get(numAcceptors*
sizeof(
Bond), (
char*)acceptors);
6180 msg->
get(numExclusions);
6183 delete [] exclusions;
6184 exclusions=
new Exclusion[numExclusions];
6185 msg->
get(numExclusions*
sizeof(
Exclusion), (
char*)exclusions);
6191 msg->
get(numConstraints);
6193 delete [] consIndexes;
6194 consIndexes =
new int32[numAtoms];
6196 msg->
get(numAtoms, consIndexes);
6200 delete [] consParams;
6201 consParams =
new ConstraintParams[numConstraints];
6203 msg->
get(numConstraints*
sizeof(ConstraintParams), (
char*)consParams);
6211 DebugM(3,
"Receiving gridforce info\n");
6213 msg->
get(numGridforceGrids);
6215 DebugM(3,
"numGridforceGrids = " << numGridforceGrids <<
"\n");
6217 delete [] numGridforces;
6218 numGridforces =
new int[numGridforceGrids];
6220 delete [] gridfrcIndexes;
6221 delete [] gridfrcParams;
6222 delete [] gridfrcGrid;
6223 gridfrcIndexes =
new int32*[numGridforceGrids];
6224 gridfrcParams =
new GridforceParams*[numGridforceGrids];
6227 int grandTotalGrids = 0;
6228 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6229 msg->
get(numGridforces[gridnum]);
6231 gridfrcIndexes[gridnum] =
new int32[numAtoms];
6232 msg->
get(numAtoms, gridfrcIndexes[gridnum]);
6234 if (numGridforces[gridnum])
6236 gridfrcParams[gridnum] =
new GridforceParams[numGridforces[gridnum]];
6237 msg->
get(numGridforces[gridnum]*
sizeof(GridforceParams), (
char*)gridfrcParams[gridnum]);
6250 msg->
get(numStirredAtoms);
6252 delete [] stirIndexes;
6253 stirIndexes =
new int32[numAtoms];
6255 msg->
get(numAtoms, stirIndexes);
6257 if (numStirredAtoms)
6259 delete [] stirParams;
6260 stirParams =
new StirParams[numStirredAtoms];
6262 msg->
get(numStirredAtoms*
sizeof(StirParams), (
char*)stirParams);
6268 msg->
get(numMovDrag);
6269 delete [] movDragIndexes;
6270 movDragIndexes =
new int32[numAtoms];
6271 msg->
get(numAtoms, movDragIndexes);
6274 delete [] movDragParams;
6275 movDragParams =
new MovDragParams[numMovDrag];
6276 msg->
get(numMovDrag*
sizeof(MovDragParams), (
char*)movDragParams);
6282 msg->
get(numRotDrag);
6283 delete [] rotDragIndexes;
6284 rotDragIndexes =
new int32[numAtoms];
6285 msg->
get(numAtoms, rotDragIndexes);
6288 delete [] rotDragParams;
6289 rotDragParams =
new RotDragParams[numRotDrag];
6290 msg->
get(numRotDrag*
sizeof(RotDragParams), (
char*)rotDragParams);
6296 msg->
get(numConsTorque);
6297 delete [] consTorqueIndexes;
6298 consTorqueIndexes =
new int32[numAtoms];
6299 msg->
get(numAtoms, consTorqueIndexes);
6302 delete [] consTorqueParams;
6303 consTorqueParams =
new ConsTorqueParams[numConsTorque];
6304 msg->
get(numConsTorque*
sizeof(ConsTorqueParams), (
char*)consTorqueParams);
6310 { msg->
get(numConsForce);
6311 delete [] consForceIndexes;
6312 consForceIndexes =
new int32[numAtoms];
6313 msg->
get(numAtoms, consForceIndexes);
6315 {
delete [] consForce;
6316 consForce =
new Vector[numConsForce];
6317 msg->
get(numConsForce*
sizeof(
Vector), (
char*)consForce);
6323 msg->
get(numMolecules);
6324 msg->
get(numLargeMolecules);
6325 delete [] moleculeAtom;
6326 delete [] moleculeStartIndex;
6327 moleculeAtom =
new int32[numAtoms];
6328 moleculeStartIndex =
new int32[numMolecules+1];
6329 msg->
get(numAtoms, moleculeAtom);
6330 msg->
get(numMolecules+1, moleculeStartIndex);
6334 exPressureAtomFlags =
new int32[numAtoms];
6335 msg->
get(numExPressureAtoms);
6336 msg->
get(numAtoms, exPressureAtomFlags);
6339 #ifndef MEM_OPT_VERSION 6343 delete [] langevinParams;
6344 langevinParams =
new Real[numAtoms];
6346 msg->
get(numAtoms, langevinParams);
6352 delete [] fixedAtomFlags;
6353 fixedAtomFlags =
new int32[numAtoms];
6355 msg->
get(numFixedAtoms);
6356 msg->
get(numAtoms, fixedAtomFlags);
6357 msg->
get(numFixedRigidBonds);
6362 if( qmAtomGroup != 0)
6363 delete [] qmAtomGroup;
6364 qmAtomGroup =
new Real[numAtoms];
6366 msg->
get(numAtoms, qmAtomGroup);
6368 msg->
get(qmNumQMAtoms);
6371 delete [] qmAtmChrg;
6372 qmAtmChrg =
new Real[qmNumQMAtoms];
6374 msg->
get(qmNumQMAtoms, qmAtmChrg);
6377 delete [] qmAtmIndx;
6378 qmAtmIndx =
new int[qmNumQMAtoms];
6380 msg->
get(qmNumQMAtoms, qmAtmIndx);
6384 msg->
get(qmNumBonds);
6386 msg->
get(qmMeNumBonds);
6388 if( qmMeMMindx != 0)
6389 delete [] qmMeMMindx;
6390 qmMeMMindx =
new int[qmMeNumBonds];
6392 msg->
get(qmMeNumBonds, qmMeMMindx);
6395 delete [] qmMeQMGrp;
6396 qmMeQMGrp =
new Real[qmMeNumBonds];
6398 msg->
get(qmMeNumBonds, qmMeQMGrp);
6402 msg->
get(qmNumGrps);
6406 qmGrpID =
new Real[qmNumGrps];
6407 msg->
get(qmNumGrps, qmGrpID);
6409 if( qmCustPCSizes != 0)
6410 delete [] qmCustPCSizes;
6411 qmCustPCSizes =
new int[qmNumGrps];
6412 msg->
get(qmNumGrps, qmCustPCSizes);
6414 msg->
get(qmTotCustPCs);
6416 if( qmCustomPCIdxs != 0)
6417 delete [] qmCustomPCIdxs;
6418 qmCustomPCIdxs =
new int[qmTotCustPCs];
6419 msg->
get(qmTotCustPCs, qmCustomPCIdxs);
6425 delete [] fepAtomFlags;
6426 fepAtomFlags =
new unsigned char[numAtoms];
6428 msg->
get(numFepInitial);
6429 msg->
get(numFepFinal);
6430 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)fepAtomFlags);
6436 delete [] ssAtomFlags;
6437 delete [] ss_vdw_type;
6439 ssAtomFlags =
new unsigned char[numAtoms];
6441 ss_index =
new int [numAtoms];
6442 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)ssAtomFlags);
6443 msg->
get(ss_num_vdw_params);
6445 msg->
get(numAtoms*
sizeof(
int), (
char*)ss_index);
6448 #ifdef OPENATOM_VERSION 6451 delete [] fepAtomFlags;
6452 fepAtomFlags =
new unsigned char[numAtoms];
6454 msg->
get(numFepInitial);
6455 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)fepAtomFlags);
6456 #endif //OPENATOM_VERSION 6459 msg->
get(is_lonepairs_psf);
6460 if (is_lonepairs_psf) {
6461 msg->
get(numLphosts);
6463 lphosts =
new Lphost[numLphosts];
6464 msg->
get(numLphosts*
sizeof(
Lphost), (
char*)lphosts);
6466 msg->
get(is_drude_psf);
6468 delete[] drudeConsts;
6471 msg->
get(numTholes);
6473 tholes =
new Thole[numTholes];
6474 msg->
get(numTholes*
sizeof(
Thole), (
char*)tholes);
6475 msg->
get(numAnisos);
6477 anisos =
new Aniso[numAnisos];
6478 msg->
get(numAnisos*
sizeof(
Aniso), (
char*)anisos);
6480 msg->
get(numZeroMassAtoms);
6485 delete [] lcpoParamType;
6486 lcpoParamType =
new int[numAtoms];
6487 msg->
get(numAtoms, (
int*)lcpoParamType);
6506 delete [] gromacsPair_type;
6509 delete [] pointerToLJBeg;
6510 pointerToLJBeg =
new int[numAtoms];
6511 msg->
get((numAtoms),pointerToLJBeg);
6512 delete [] pointerToLJEnd;
6513 pointerToLJEnd =
new int[numAtoms];
6514 msg->
get((numAtoms),pointerToLJEnd);
6527 delete [] indxGaussA;
6530 delete [] indxGaussB;
6548 delete [] gRepulsive;
6551 delete [] pointerToGaussBeg;
6552 pointerToGaussBeg =
new int[numAtoms];
6553 msg->
get((numAtoms),pointerToGaussBeg);
6554 delete [] pointerToGaussEnd;
6555 pointerToGaussEnd =
new int[numAtoms];
6556 msg->
get((numAtoms),pointerToGaussEnd);
6564 #ifdef MEM_OPT_VERSION 6566 build_excl_check_signatures();
6569 numBonds = numCalcBonds = 0;
6570 numAngles = numCalcAngles = 0;
6571 numDihedrals = numCalcDihedrals = 0;
6572 numImpropers = numCalcImpropers = 0;
6573 numCrossterms = numCalcCrossterms = 0;
6574 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;
6582 build_atom_status();
6583 build_lists_by_atom();
6622 DebugM(3,
"Entered build_gridforce_params multi...\n");
6627 numGridforceGrids = 0;
6628 while (mgridParams != NULL) {
6629 numGridforceGrids++;
6630 mgridParams = mgridParams->
next;
6633 DebugM(3,
"numGridforceGrids = " << numGridforceGrids <<
"\n");
6634 gridfrcIndexes =
new int32*[numGridforceGrids];
6635 gridfrcParams =
new GridforceParams*[numGridforceGrids];
6637 numGridforces =
new int[numGridforceGrids];
6639 int grandTotalGrids = 0;
6641 mgridParams =
simParams->mgridforcelist.get_first();
6642 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6643 int current_index=0;
6650 if (mgridParams == NULL) {
6651 NAMD_die(
"Problem with mgridParams!");
6657 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, gridforceFile required.");
6664 if ( (cwd == NULL) || (mgridParams->
gridforceFile[0] ==
'/') )
6670 strcpy(filename, cwd);
6674 kPDB =
new PDB(filename);
6677 NAMD_die(
"Memory allocation failed in Molecule::build_gridforce_params");
6682 NAMD_die(
"Number of atoms in grid force PDB doesn't match coordinate PDB");
6701 else if (strcasecmp(mgridParams->
gridforceCol,
"Y") == 0)
6705 else if (strcasecmp(mgridParams->
gridforceCol,
"Z") == 0)
6709 else if (strcasecmp(mgridParams->
gridforceCol,
"O") == 0)
6713 else if (strcasecmp(mgridParams->
gridforceCol,
"B") == 0)
6719 NAMD_die(
"gridforcecol must have value of X, Y, Z, O, or B");
6752 NAMD_die(
"gridforcechargecol must have value of X, Y, Z, O, or B");
6757 NAMD_die(
"gridforcecol and gridforcechargecol cannot have same value");
6764 gridfrcIndexes[gridnum] =
new int32[numAtoms];
6766 if (gridfrcIndexes[gridnum] == NULL)
6768 NAMD_die(
"memory allocation failed in Molecule::build_gridforce_params()");
6772 for (i=0; i<numAtoms; i++)
6778 kval = (kPDB->
atom(i))->xcoor();
6781 kval = (kPDB->
atom(i))->ycoor();
6784 kval = (kPDB->
atom(i))->zcoor();
6787 kval = (kPDB->
atom(i))->occupancy();
6790 kval = (kPDB->
atom(i))->temperaturefactor();
6797 gridfrcIndexes[gridnum][i] = current_index;
6803 gridfrcIndexes[gridnum][i] = -1;
6807 if (current_index == 0)
6810 iout <<
iWARN <<
"NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" <<
endi;
6815 gridfrcParams[gridnum] =
new GridforceParams[current_index];
6816 if (gridfrcParams[gridnum] == NULL)
6818 NAMD_die(
"memory allocation failed in Molecule::build_gridforce_params");
6822 numGridforces[gridnum] = current_index;
6826 for (i=0; i<numAtoms; i++)
6828 if (gridfrcIndexes[gridnum][i] != -1)
6834 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->xcoor();
6837 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->ycoor();
6840 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->zcoor();
6843 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->occupancy();
6846 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->temperaturefactor();
6854 #ifdef MEM_OPT_VERSION 6855 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
6857 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
6861 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->xcoor();
6864 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->ycoor();
6867 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->zcoor();
6870 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->occupancy();
6873 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->temperaturefactor();
6894 strcpy(potfilename, cwd);
6901 DebugM(3,
"allocating GridforceGrid(" << gridnum <<
")\n");
6905 DebugM(4,
"grandTotalGrids = " << grandTotalGrids <<
"\n" <<
endi);
6908 mgridParams = mgridParams->
next;
6913 #ifdef MEM_OPT_VERSION 6914 void Molecule::delEachAtomSigs(){
6920 if(CmiMyRank())
return;
6922 delete [] eachAtomSig;
6923 delete [] eachAtomExclSig;
6925 eachAtomExclSig = NULL;
6928 void Molecule::delChargeSpace(){
6929 if(CmiMyRank())
return;
6931 delete [] atomChargePool;
6932 delete [] eachAtomCharge;
6933 atomChargePool = NULL;
6934 eachAtomCharge = NULL;
6937 void Molecule::delMassSpace(){
6938 if(CmiMyRank())
return;
6940 delete [] atomMassPool;
6941 delete [] eachAtomMass;
6942 atomMassPool = NULL;
6943 eachAtomMass = NULL;
6946 void Molecule::delClusterSigs() {
6947 if(CmiMyRank())
return;
6949 delete [] clusterSigs;
6953 void Molecule::delAtomNames(){
6954 if(CmiMyRank())
return;
6956 delete [] atomNames;
6961 void Molecule::delFixedAtoms(){
6962 if(CmiMyRank())
return;
6963 delete fixedAtomsSet;
6964 fixedAtomsSet = NULL;
6969 #endif // MOLECULE2_C undefined = first object file 6970 #ifdef MOLECULE2_C // second object file 7002 int current_index=0;
7012 if (consref == NULL)
7014 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, consref required.");
7015 refPDB = initial_pdb;
7019 if (consref->
next != NULL)
7021 NAMD_die(
"Multiple definitions of constraint reference file in configruation file");
7024 if ( (cwd == NULL) || (consref->
data[0] ==
'/') )
7026 strcpy(filename, consref->
data);
7030 strcpy(filename, cwd);
7031 strcat(filename, consref->
data);
7034 refPDB =
new PDB(filename);
7035 if ( refPDB == NULL )
7037 NAMD_die(
"Memory allocation failed in Molecule::build_constraint_params");
7042 NAMD_die(
"Number of atoms in constraint reference PDB doesn't match coordinate PDB");
7049 if (conskfile == NULL)
7051 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, conskfile required.");
7056 if (conskfile->
next != NULL)
7058 NAMD_die(
"Multiple definitions of constraint constant file in configuration file");
7061 if ( (consref != NULL) && (strcasecmp(consref->
data, conskfile->
data) == 0) )
7068 if ( (cwd == NULL) || (conskfile->
data[0] ==
'/') )
7070 strcpy(filename, conskfile->
data);
7074 strcpy(filename, cwd);
7075 strcat(filename, conskfile->
data);
7078 kPDB =
new PDB(filename);
7081 NAMD_die(
"Memory allocation failed in Molecule::build_constraint_params");
7086 NAMD_die(
"Number of atoms in constraint constant PDB doesn't match coordinate PDB");
7096 if (conskcol == NULL)
7102 if (conskcol->
next != NULL)
7104 NAMD_die(
"Multiple definitions of harmonic constraint column in config file");
7107 if (strcasecmp(conskcol->
data,
"X") == 0)
7111 else if (strcasecmp(conskcol->
data,
"Y") == 0)
7115 else if (strcasecmp(conskcol->
data,
"Z") == 0)
7119 else if (strcasecmp(conskcol->
data,
"O") == 0)
7123 else if (strcasecmp(conskcol->
data,
"B") == 0)
7129 NAMD_die(
"conskcol must have value of X, Y, Z, O, or B");
7136 consIndexes =
new int32[numAtoms];
7138 if (consIndexes == NULL)
7140 NAMD_die(
"memory allocation failed in Molecule::build_constraint_params()");
7144 for (i=0; i<numAtoms; i++)
7150 kval = (kPDB->
atom(i))->xcoor();
7153 kval = (kPDB->
atom(i))->ycoor();
7156 kval = (kPDB->
atom(i))->zcoor();
7159 kval = (kPDB->
atom(i))->occupancy();
7162 kval = (kPDB->
atom(i))->temperaturefactor();
7169 consIndexes[i] = current_index;
7175 consIndexes[i] = -1;
7179 if (current_index == 0)
7182 iout <<
iWARN <<
"NO CONSTRAINED ATOMS WERE FOUND, BUT CONSTRAINTS ARE ON . . .\n" <<
endi;
7187 consParams =
new ConstraintParams[current_index];
7189 if (consParams == NULL)
7191 NAMD_die(
"memory allocation failed in Molecule::build_constraint_params");
7195 numConstraints = current_index;
7199 for (i=0; i<numAtoms; i++)
7201 if (consIndexes[i] != -1)
7207 consParams[consIndexes[i]].k = (kPDB->
atom(i))->xcoor();
7210 consParams[consIndexes[i]].k = (kPDB->
atom(i))->ycoor();
7213 consParams[consIndexes[i]].k = (kPDB->
atom(i))->zcoor();
7216 consParams[consIndexes[i]].k = (kPDB->
atom(i))->occupancy();
7219 consParams[consIndexes[i]].k = (kPDB->
atom(i))->temperaturefactor();
7224 consParams[consIndexes[i]].refPos.x = (refPDB->
atom(i))->xcoor();
7225 consParams[consIndexes[i]].refPos.y = (refPDB->
atom(i))->ycoor();
7226 consParams[consIndexes[i]].refPos.z = (refPDB->
atom(i))->zcoor();
7231 if (consref != NULL)
7236 if ((conskfile != NULL) &&
7237 !((consref != NULL) &&
7238 (strcasecmp(consref->
data, conskfile->
data) == 0)
7277 int current_index=0;
7286 if (movDragFile == NULL) {
7287 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, movDragFile required.");
7292 if (movDragFile->
next != NULL) {
7293 NAMD_die(
"Multiple definitions of moving drag tag file in configuration file");
7296 if ( (cwd == NULL) || (movDragFile->
data[0] ==
'/') ) {
7297 strcpy(mainfilename, movDragFile->
data);
7299 strcpy(mainfilename, cwd);
7300 strcat(mainfilename, movDragFile->
data);
7303 tPDB =
new PDB(mainfilename);
7304 if ( tPDB == NULL ) {
7305 NAMD_die(
"Memory allocation failed in Molecule::build_movdrag_params");
7309 NAMD_die(
"Number of atoms in moving drag tag PDB doesn't match coordinate PDB");
7317 if (movDragVelFile == NULL) {
7318 if (movDragFile == NULL) {
7319 NAMD_die(
"Moving drag velocity file can not be same as coordinate PDB file");
7321 if (movDragVelFile->
next != NULL) {
7322 NAMD_die(
"Multiple definitions of moving drag velocity file in configuration file");
7329 if ( (cwd == NULL) || (movDragVelFile->
data[0] ==
'/') ) {
7330 strcpy(velfilename, movDragVelFile->
data);
7332 strcpy(velfilename, cwd);
7333 strcat(velfilename, movDragVelFile->
data);
7336 vPDB =
new PDB(velfilename);
7337 if ( vPDB == NULL ) {
7338 NAMD_die(
"Memory allocation failed in Molecule::build_movdrag_params");
7342 NAMD_die(
"Number of atoms in moving drag velocity PDB doesn't match coordinate PDB");
7354 if (movDragCol == NULL) {
7357 if (movDragCol->
next != NULL) {
7358 NAMD_die(
"Multiple definitions of drag column in config file");
7361 if (movDragFile == NULL
7362 && strcasecmp(movDragCol->
data,
"B")
7363 && strcasecmp(movDragCol->
data,
"O")) {
7364 NAMD_die(
"Can not read moving drag tags from X, Y, or Z column of the coordinate or velocity file");
7366 if (!strcasecmp(movDragCol->
data,
"X")) {
7368 }
else if (!strcasecmp(movDragCol->
data,
"Y")) {
7370 }
else if (!strcasecmp(movDragCol->
data,
"Z")) {
7372 }
else if (!strcasecmp(movDragCol->
data,
"O")) {
7374 }
else if (!strcasecmp(movDragCol->
data,
"B")) {
7378 NAMD_die(
"movDragCol must have value of X, Y, Z, O, or B");
7385 movDragIndexes =
new int32[numAtoms];
7386 if (movDragIndexes == NULL) {
7387 NAMD_die(
"memory allocation failed in Molecule::build_movdrag_params()");
7391 for (i=0; i<numAtoms; i++) {
7394 dtval = (tPDB->
atom(i))->xcoor();
7397 dtval = (tPDB->
atom(i))->ycoor();
7400 dtval = (tPDB->
atom(i))->zcoor();
7403 dtval = (tPDB->
atom(i))->occupancy();
7406 dtval = (tPDB->
atom(i))->temperaturefactor();
7412 movDragIndexes[i] = current_index;
7416 movDragIndexes[i] = -1;
7420 if (current_index == 0) {
7422 iout <<
iWARN <<
"NO DRAGGED ATOMS WERE FOUND, BUT MOVING DRAG IS ON . . . " <<
endi;
7425 movDragParams =
new MovDragParams[current_index];
7426 if (movDragParams == NULL) {
7427 NAMD_die(
"memory allocation failed in Molecule::build_movdrag_params");
7431 numMovDrag = current_index;
7435 for (i=0; i<numAtoms; i++) {
7436 if (movDragIndexes[i] != -1) {
7437 movDragParams[movDragIndexes[i]].v[0] = (vPDB->
atom(i))->xcoor();
7438 movDragParams[movDragIndexes[i]].v[1] = (vPDB->
atom(i))->ycoor();
7439 movDragParams[movDragIndexes[i]].v[2] = (vPDB->
atom(i))->zcoor();
7443 if (movDragFile != NULL)
delete tPDB;
7444 if (movDragVelFile != NULL)
delete vPDB;
7481 PDB *tPDB, *aPDB, *pPDB, *vPDB;
7483 int current_index=0;
7496 if (rotDragFile == NULL) {
7497 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, rotDragFile required.");
7502 if (rotDragFile->
next != NULL) {
7503 NAMD_die(
"Multiple definitions of rotating drag tag file in configuration file");
7506 if ( (cwd == NULL) || (rotDragFile->
data[0] ==
'/') ) {
7507 strcpy(mainfilename, rotDragFile->
data);
7509 strcpy(mainfilename, cwd);
7510 strcat(mainfilename, rotDragFile->
data);
7513 tPDB =
new PDB(mainfilename);
7514 if ( tPDB == NULL ) {
7515 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7519 NAMD_die(
"Number of atoms in rotating drag tag PDB doesn't match coordinate PDB");
7527 if (rotDragAxisFile == NULL) {
7528 if (rotDragFile == NULL) {
7529 NAMD_die(
"Rotating drag axis file can not be same as coordinate PDB file");
7531 if (rotDragAxisFile->
next != NULL) {
7532 NAMD_die(
"Multiple definitions of rotating drag axis file in configuration file");
7534 if (rotDragPivotFile == NULL) {
7535 NAMD_die(
"Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
7542 if ( (cwd == NULL) || (rotDragAxisFile->
data[0] ==
'/') ) {
7543 strcpy(axisfilename, rotDragAxisFile->
data);
7545 strcpy(axisfilename, cwd);
7546 strcat(axisfilename, rotDragAxisFile->
data);
7549 aPDB =
new PDB(axisfilename);
7550 if ( aPDB == NULL ) {
7551 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7555 NAMD_die(
"Number of atoms in rotating drag axis PDB doesn't match coordinate PDB");
7563 if (rotDragPivotFile == NULL) {
7564 if (rotDragFile == NULL) {
7565 NAMD_die(
"Rotating drag pivot point file can not be same as coordinate PDB file");
7567 if (rotDragPivotFile->
next != NULL) {
7568 NAMD_die(
"Multiple definitions of rotating drag pivot point file in configuration file");
7570 if (rotDragAxisFile == NULL) {
7571 NAMD_die(
"Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
7578 if ( (cwd == NULL) || (rotDragPivotFile->
data[0] ==
'/') ) {
7579 strcpy(pivotfilename, rotDragPivotFile->
data);
7581 strcpy(pivotfilename, cwd);
7582 strcat(pivotfilename, rotDragPivotFile->
data);
7585 pPDB =
new PDB(pivotfilename);
7586 if ( pPDB == NULL ) {
7587 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7591 NAMD_die(
"Number of atoms in rotating drag pivot point PDB doesn't match coordinate PDB");
7600 if (rotDragVelFile == NULL) {
7603 if (rotDragVelFile->
next != NULL) {
7604 NAMD_die(
"Multiple definitions of rotating drag velocity file in configuration file");
7607 if ( (cwd == NULL) || (rotDragVelFile->
data[0] ==
'/') ) {
7608 strcpy(velfilename, rotDragVelFile->
data);
7610 strcpy(velfilename, cwd);
7611 strcat(velfilename, rotDragVelFile->
data);
7614 vPDB =
new PDB(velfilename);
7615 if ( vPDB == NULL ) {
7616 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7620 NAMD_die(
"Number of atoms in rotating drag velocity PDB doesn't match coordinate PDB");
7631 if (rotDragCol == NULL) {
7634 if (rotDragCol->
next != NULL) {
7635 NAMD_die(
"Multiple definitions of drag tag column in config file");
7638 if ( rotDragFile == NULL
7639 && (!strcasecmp(rotDragCol->
data,
"X")
7640 || !strcasecmp(rotDragCol->
data,
"Y")
7641 || !strcasecmp(rotDragCol->
data,
"Z"))) {
7642 NAMD_die(
"Can not read rotating drag tags from X, Y, or Z column of the PDB coordinate file");
7644 if (!strcasecmp(rotDragCol->
data,
"X")) {
7646 }
else if (!strcasecmp(rotDragCol->
data,
"Y")) {
7648 }
else if (!strcasecmp(rotDragCol->
data,
"Z")) {
7650 }
else if (!strcasecmp(rotDragCol->
data,
"O")) {
7652 }
else if (!strcasecmp(rotDragCol->
data,
"B")) {
7656 NAMD_die(
"rotDragCol must have value of X, Y, Z, O, or B");
7668 if (rotDragVelCol == NULL) {
7671 if (rotDragVelCol->
next != NULL) {
7672 NAMD_die(
"Multiple definitions of drag angular velocity column in config file");
7675 if (rotDragVelFile == NULL
7676 && rotDragFile == NULL
7677 && strcasecmp(rotDragCol->
data,
"B")
7678 && strcasecmp(rotDragCol->
data,
"O")) {
7679 NAMD_die(
"Can not read rotating drag angular velocities from X, Y, or Z column of the PDB coordinate file");
7681 if (!strcasecmp(rotDragVelCol->
data,
"X")) {
7683 }
else if (!strcasecmp(rotDragVelCol->
data,
"Y")) {
7685 }
else if (!strcasecmp(rotDragVelCol->
data,
"Z")) {
7687 }
else if (!strcasecmp(rotDragVelCol->
data,
"O")) {
7689 }
else if (!strcasecmp(rotDragVelCol->
data,
"B")) {
7693 NAMD_die(
"rotDragVelCol must have value of X, Y, Z, O, or B");
7700 rotDragIndexes =
new int32[numAtoms];
7701 if (rotDragIndexes == NULL) {
7702 NAMD_die(
"memory allocation failed in Molecule::build_rotdrag_params()");
7706 for (i=0; i<numAtoms; i++) {
7709 dtval = (tPDB->
atom(i))->xcoor();
7712 dtval = (tPDB->
atom(i))->ycoor();
7715 dtval = (tPDB->
atom(i))->zcoor();
7718 dtval = (tPDB->
atom(i))->occupancy();
7721 dtval = (tPDB->
atom(i))->temperaturefactor();
7727 rotDragIndexes[i] = current_index;
7731 rotDragIndexes[i] = -1;
7735 if (current_index == 0) {
7736 iout <<
iWARN <<
"NO DRAGGED ATOMS WERE FOUND, BUT ROTATING DRAG IS ON . . . " <<
endi;
7738 rotDragParams =
new RotDragParams[current_index];
7739 if (rotDragParams == NULL) {
7740 NAMD_die(
"memory allocation failed in Molecule::build_rotdrag_params");
7744 numRotDrag = current_index;
7748 for (i=0; i<numAtoms; i++) {
7749 if (rotDragIndexes[i] != -1) {
7750 rotDragParams[rotDragIndexes[i]].a[0] = (aPDB->
atom(i))->xcoor();
7751 rotDragParams[rotDragIndexes[i]].a[1] = (aPDB->
atom(i))->ycoor();
7752 rotDragParams[rotDragIndexes[i]].a[2] = (aPDB->
atom(i))->zcoor();
7753 rotDragParams[rotDragIndexes[i]].p[0] = (pPDB->
atom(i))->xcoor();
7754 rotDragParams[rotDragIndexes[i]].p[1] = (pPDB->
atom(i))->ycoor();
7755 rotDragParams[rotDragIndexes[i]].p[2] = (pPDB->
atom(i))->zcoor();
7758 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->xcoor();
7761 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->ycoor();
7764 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->zcoor();
7767 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->occupancy();
7770 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->temperaturefactor();
7776 if (rotDragFile != NULL)
delete tPDB;
7777 if (rotDragAxisFile != NULL)
delete aPDB;
7778 if (rotDragPivotFile != NULL)
delete pPDB;
7779 if (rotDragVelFile != NULL)
delete vPDB;
7816 PDB *tPDB, *aPDB, *pPDB, *vPDB;
7818 int current_index=0;
7831 if (consTorqueFile == NULL) {
7832 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, consTorqueFile required.");
7837 if (consTorqueFile->
next != NULL) {
7838 NAMD_die(
"Multiple definitions of \"constant\" torque tag file in configuration file");
7841 if ( (cwd == NULL) || (consTorqueFile->
data[0] ==
'/') ) {
7842 strcpy(mainfilename, consTorqueFile->
data);
7844 strcpy(mainfilename, cwd);
7845 strcat(mainfilename, consTorqueFile->
data);
7848 tPDB =
new PDB(mainfilename);
7849 if ( tPDB == NULL ) {
7850 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7854 NAMD_die(
"Number of atoms in \"constant\" torque tag PDB doesn't match coordinate PDB");
7862 if (consTorqueAxisFile == NULL) {
7863 if (consTorqueFile == NULL) {
7864 NAMD_die(
"\"Constant\" torque axis file can not be same as coordinate PDB file");
7866 if (consTorqueAxisFile->
next != NULL) {
7867 NAMD_die(
"Multiple definitions of \"constant\" torque axis file in configuration file");
7869 if (consTorquePivotFile == NULL) {
7870 NAMD_die(
"Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
7877 if ( (cwd == NULL) || (consTorqueAxisFile->
data[0] ==
'/') ) {
7878 strcpy(axisfilename, consTorqueAxisFile->
data);
7880 strcpy(axisfilename, cwd);
7881 strcat(axisfilename, consTorqueAxisFile->
data);
7884 aPDB =
new PDB(axisfilename);
7885 if ( aPDB == NULL ) {
7886 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7890 NAMD_die(
"Number of atoms in \"constant\" torque axis PDB doesn't match coordinate PDB");
7898 if (consTorquePivotFile == NULL) {
7899 if (consTorqueFile == NULL) {
7900 NAMD_die(
"\"Constant\" torque pivot point file can not be same as coordinate PDB file");
7902 if (consTorquePivotFile->
next != NULL) {
7903 NAMD_die(
"Multiple definitions of \"constant\" torque pivot point file in configuration file");
7905 if (consTorqueAxisFile == NULL) {
7906 NAMD_die(
"Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
7913 if ( (cwd == NULL) || (consTorquePivotFile->
data[0] ==
'/') ) {
7914 strcpy(pivotfilename, consTorquePivotFile->
data);
7916 strcpy(pivotfilename, cwd);
7917 strcat(pivotfilename, consTorquePivotFile->
data);
7920 pPDB =
new PDB(pivotfilename);
7921 if ( pPDB == NULL ) {
7922 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7926 NAMD_die(
"Number of atoms in \"constant\" torque pivot point PDB doesn't match coordinate PDB");
7935 if (consTorqueValFile == NULL) {
7938 if (consTorqueValFile->
next != NULL) {
7939 NAMD_die(
"Multiple definitions of \"constant\" torque velocity file in configuration file");
7942 if ( (cwd == NULL) || (consTorqueValFile->
data[0] ==
'/') ) {
7943 strcpy(velfilename, consTorqueValFile->
data);
7945 strcpy(velfilename, cwd);
7946 strcat(velfilename, consTorqueValFile->
data);
7949 vPDB =
new PDB(velfilename);
7950 if ( vPDB == NULL ) {
7951 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7955 NAMD_die(
"Number of atoms in \"constant\" torque velocity PDB doesn't match coordinate PDB");
7966 if (consTorqueCol == NULL) {
7969 if (consTorqueCol->
next != NULL) {
7970 NAMD_die(
"Multiple definitions of torque tag column in config file");
7973 if ( consTorqueFile == NULL
7974 && (!strcasecmp(consTorqueCol->
data,
"X")
7975 || !strcasecmp(consTorqueCol->
data,
"Y")
7976 || !strcasecmp(consTorqueCol->
data,
"Z"))) {
7977 NAMD_die(
"Can not read \"constant\" torque tags from X, Y, or Z column of the PDB coordinate file");
7979 if (!strcasecmp(consTorqueCol->
data,
"X")) {
7981 }
else if (!strcasecmp(consTorqueCol->
data,
"Y")) {
7983 }
else if (!strcasecmp(consTorqueCol->
data,
"Z")) {
7985 }
else if (!strcasecmp(consTorqueCol->
data,
"O")) {
7987 }
else if (!strcasecmp(consTorqueCol->
data,
"B")) {
7991 NAMD_die(
"consTorqueCol must have value of X, Y, Z, O, or B");
8003 if (consTorqueValCol == NULL) {
8006 if (consTorqueValCol->
next != NULL) {
8007 NAMD_die(
"Multiple definitions of torque value column in config file");
8010 if (consTorqueValFile == NULL
8011 && consTorqueFile == NULL
8012 && strcasecmp(consTorqueCol->
data,
"B")
8013 && strcasecmp(consTorqueCol->
data,
"O")) {
8014 NAMD_die(
"Can not read \"constant\" torque values from X, Y, or Z column of the PDB coordinate file");
8016 if (!strcasecmp(consTorqueValCol->
data,
"X")) {
8018 }
else if (!strcasecmp(consTorqueValCol->
data,
"Y")) {
8020 }
else if (!strcasecmp(consTorqueValCol->
data,
"Z")) {
8022 }
else if (!strcasecmp(consTorqueValCol->
data,
"O")) {
8024 }
else if (!strcasecmp(consTorqueValCol->
data,
"B")) {
8028 NAMD_die(
"consTorqueValCol must have value of X, Y, Z, O, or B");
8035 consTorqueIndexes =
new int32[numAtoms];
8036 if (consTorqueIndexes == NULL) {
8037 NAMD_die(
"memory allocation failed in Molecule::build_constorque_params()");
8041 for (i=0; i<numAtoms; i++) {
8044 dtval = (tPDB->
atom(i))->xcoor();
8047 dtval = (tPDB->
atom(i))->ycoor();
8050 dtval = (tPDB->
atom(i))->zcoor();
8053 dtval = (tPDB->
atom(i))->occupancy();
8056 dtval = (tPDB->
atom(i))->temperaturefactor();
8062 consTorqueIndexes[i] = current_index;
8066 consTorqueIndexes[i] = -1;
8070 if (current_index == 0) {
8071 iout <<
iWARN <<
"NO TORQUED ATOMS WERE FOUND, BUT \"CONSTANT\" TORQUE IS ON . . . " <<
endi;
8073 consTorqueParams =
new ConsTorqueParams[current_index];
8074 if (consTorqueParams == NULL) {
8075 NAMD_die(
"memory allocation failed in Molecule::build_constorque_params");
8079 numConsTorque = current_index;
8083 for (i=0; i<numAtoms; i++) {
8084 if (consTorqueIndexes[i] != -1) {
8085 consTorqueParams[consTorqueIndexes[i]].a[0] = (aPDB->
atom(i))->xcoor();
8086 consTorqueParams[consTorqueIndexes[i]].a[1] = (aPDB->
atom(i))->ycoor();
8087 consTorqueParams[consTorqueIndexes[i]].a[2] = (aPDB->
atom(i))->zcoor();
8088 consTorqueParams[consTorqueIndexes[i]].p[0] = (pPDB->
atom(i))->xcoor();
8089 consTorqueParams[consTorqueIndexes[i]].p[1] = (pPDB->
atom(i))->ycoor();
8090 consTorqueParams[consTorqueIndexes[i]].p[2] = (pPDB->
atom(i))->zcoor();
8093 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->xcoor();
8096 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->ycoor();
8099 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->zcoor();
8102 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->occupancy();
8105 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->temperaturefactor();
8111 if (consTorqueFile != NULL)
delete tPDB;
8112 if (consTorqueAxisFile != NULL)
delete aPDB;
8113 if (consTorquePivotFile != NULL)
delete pPDB;
8114 if (consTorqueValFile != NULL)
delete vPDB;
8140 iout <<
iWARN <<
"NO CONSTANT FORCES SPECIFIED, BUT CONSTANT FORCE IS ON . . .\n" <<
endi;
8141 consForceIndexes =
new int32[numAtoms];
8142 for (i=0; i<numAtoms; i++) consForceIndexes[i] = -1;
8146 if ((forcePDB=
new PDB(filename)) == NULL)
8147 NAMD_die(
"Memory allocation failed in Molecule::build_constant_forces");
8149 NAMD_die(
"Number of atoms in constant force PDB doesn't match coordinate PDB");
8154 consForceIndexes =
new int32[numAtoms];
8155 if (consForceIndexes == NULL)
8156 NAMD_die(
"memory allocation failed in Molecule::build_constant_forces()");
8160 for (i=0; i<numAtoms; i++)
8164 consForceIndexes[i] = -1;
8167 consForceIndexes[i] = numConsForce++;
8169 if (numConsForce == 0)
8171 iout <<
iWARN <<
"NO NON-ZERO FORCES WERE FOUND, BUT CONSTANT FORCE IS ON . . .\n" <<
endi;
8174 consForce =
new Vector[numConsForce];
8175 if (consForce == NULL)
8176 NAMD_die(
"memory allocation failed in Molecule::build_constant_forces");
8178 for (i=0; i<numAtoms; i++)
8179 if ((index=consForceIndexes[i]) != -1)
8196 langevinParams =
new Real[numAtoms];
8198 if ( (langevinParams == NULL) )
8200 NAMD_die(
"memory allocation failed in Molecule::build_langevin_params()");
8204 for (
int i=0; i<numAtoms; i++)
8208 if ( (! doHydrogen) && is_hydrogen(i) ) bval = 0;
8209 else if ( is_lp(i) ) bval = 0;
8210 else if ( is_drude(i) ) bval = drudeCoupling;
8213 langevinParams[i] = bval;
8250 if (langfile == NULL)
8252 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, langevinFile required.");
8257 if (langfile->
next != NULL)
8259 NAMD_die(
"Multiple definitions of langvein PDB file in configuration file");
8262 if ( (cwd == NULL) || (langfile->
data[0] ==
'/') )
8264 strcpy(filename, langfile->
data);
8268 strcpy(filename, cwd);
8269 strcat(filename, langfile->
data);
8272 bPDB =
new PDB(filename);
8275 NAMD_die(
"Memory allocation failed in Molecule::build_langevin_params");
8280 NAMD_die(
"Number of atoms in langevin parameter PDB doesn't match coordinate PDB");
8289 if (langcol == NULL)
8295 if (langcol->
next != NULL)
8297 NAMD_die(
"Multiple definitions of langevin parameter column in config file");
8300 if (strcasecmp(langcol->
data,
"X") == 0)
8304 else if (strcasecmp(langcol->
data,
"Y") == 0)
8308 else if (strcasecmp(langcol->
data,
"Z") == 0)
8312 else if (strcasecmp(langcol->
data,
"O") == 0)
8316 else if (strcasecmp(langcol->
data,
"B") == 0)
8322 NAMD_die(
"langevincol must have value of X, Y, Z, O, or B");
8327 langevinParams =
new Real[numAtoms];
8329 if ( (langevinParams == NULL) )
8331 NAMD_die(
"memory allocation failed in Molecule::build_langevin_params()");
8335 for (i=0; i<numAtoms; i++)
8341 bval = (bPDB->
atom(i))->xcoor();
8344 bval = (bPDB->
atom(i))->ycoor();
8347 bval = (bPDB->
atom(i))->zcoor();
8350 bval = (bPDB->
atom(i))->occupancy();
8353 bval = (bPDB->
atom(i))->temperaturefactor();
8358 langevinParams[i] = bval;
8362 if (langfile != NULL)
8401 if (fixedfile == NULL)
8403 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, fixedAtomsFile required.");
8408 if (fixedfile->
next != NULL)
8410 NAMD_die(
"Multiple definitions of fixed atoms PDB file in configuration file");
8413 if ( (cwd == NULL) || (fixedfile->
data[0] ==
'/') )
8415 strcpy(filename, fixedfile->
data);
8419 strcpy(filename, cwd);
8420 strcat(filename, fixedfile->
data);
8423 bPDB =
new PDB(filename);
8426 NAMD_die(
"Memory allocation failed in Molecule::build_fixed_atoms");
8431 NAMD_die(
"Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
8440 if (fixedcol == NULL)
8446 if (fixedcol->
next != NULL)
8448 NAMD_die(
"Multiple definitions of fixed atoms column in config file");
8451 if (strcasecmp(fixedcol->
data,
"X") == 0)
8455 else if (strcasecmp(fixedcol->
data,
"Y") == 0)
8459 else if (strcasecmp(fixedcol->
data,
"Z") == 0)
8463 else if (strcasecmp(fixedcol->
data,
"O") == 0)
8467 else if (strcasecmp(fixedcol->
data,
"B") == 0)
8473 NAMD_die(
"fixedatomscol must have value of X, Y, Z, O, or B");
8478 fixedAtomFlags =
new int32[numAtoms];
8480 if (fixedAtomFlags == NULL)
8482 NAMD_die(
"memory allocation failed in Molecule::build_fixed_atoms()");
8488 for (i=0; i<numAtoms; i++)
8494 bval = (bPDB->
atom(i))->xcoor();
8497 bval = (bPDB->
atom(i))->ycoor();
8500 bval = (bPDB->
atom(i))->zcoor();
8503 bval = (bPDB->
atom(i))->occupancy();
8506 bval = (bPDB->
atom(i))->temperaturefactor();
8512 fixedAtomFlags[i] = 1;
8516 fixedAtomFlags[i] = 0;
8521 if (fixedfile != NULL)
8528 if ( numRigidBonds ) {
8530 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
8531 int parentIsFixed = 0;
8532 for( ; h_i != h_e; ++h_i ) {
8534 parentIsFixed = fixedAtomFlags[h_i->
atomID];
8535 if ( (rigidBondLengths[h_i->
atomID] > 0.)
8536 && fixedAtomFlags[h_i[1].atomID]
8537 && fixedAtomFlags[h_i[2].atomID] ) {
8538 ++numFixedRigidBonds;
8541 if ( (rigidBondLengths[h_i->
atomID] > 0.)
8542 && fixedAtomFlags[h_i->
atomID]
8543 && parentIsFixed ) {
8544 ++numFixedRigidBonds;
8554 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
8556 for( ; h_i != h_e; ++h_i ) {
8558 if ( allFixed ) ++numFixedGroups;
8561 allFixed = allFixed && fixedAtomFlags[h_i->
atomID];
8563 if ( allFixed ) ++numFixedGroups;
8600 int current_index=0;
8608 if (stirredfile == NULL)
8610 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, stirFilename required.");
8613 iout <<
iWARN <<
"STIRRING USING INITIAL POSITION FILE FOR REFERENCE POSITIONS" <<
endi;
8617 if (stirredfile->
next != NULL)
8619 NAMD_die(
"Multiple definitions of stirred atoms PDB file in configuration file");
8622 if ( (cwd == NULL) || (stirredfile->
data[0] ==
'/') )
8624 strcpy(filename, stirredfile->
data);
8628 strcpy(filename, cwd);
8629 strcat(filename, stirredfile->
data);
8633 sPDB =
new PDB(filename);
8637 NAMD_die(
"Memory allocation failed in Molecule::build_stirred_atoms");
8643 NAMD_die(
"Number of atoms in stirred atoms PDB doesn't match coordinate PDB");
8655 if (stirredcol == NULL)
8661 if (stirredcol->
next != NULL)
8663 NAMD_die(
"Multiple definitions of stirred atoms column in config file");
8666 if (strcasecmp(stirredcol->
data,
"O") == 0)
8670 else if (strcasecmp(stirredcol->
data,
"B") == 0)
8676 NAMD_die(
"stirredAtomsCol must have value of O or B");
8683 stirIndexes =
new int32[numAtoms];
8685 if (stirIndexes == NULL)
8687 NAMD_die(
"memory allocation failed in Molecule::build_stirred_params()");
8692 for (i=0; i<numAtoms; i++)
8702 bval = (sPDB->
atom(i))->occupancy();
8705 bval = (sPDB->
atom(i))->temperaturefactor();
8714 stirIndexes[i] = current_index;
8720 stirIndexes[i] = -1;
8728 if (current_index == 0)
8731 iout <<
iWARN <<
"NO STIRRED ATOMS WERE FOUND, BUT STIRRING TORQUES ARE ON . . .\n" <<
endi;
8736 stirParams =
new StirParams[current_index];
8738 if (stirParams == NULL)
8740 NAMD_die(
"memory allocation failed in Molecule::build_stir_params");
8744 numStirredAtoms = current_index;
8748 for (i=0; i<numAtoms; i++)
8750 if (stirIndexes[i] != -1)
8754 stirParams[stirIndexes[i]].refPos.x = (sPDB->
atom(i))->xcoor();
8755 stirParams[stirIndexes[i]].refPos.y = (sPDB->
atom(i))->ycoor();
8756 stirParams[stirIndexes[i]].refPos.z = (sPDB->
atom(i))->zcoor();
8761 if (stirredfile != NULL)
8777 int a1,a2,a3,a4;
float k, ref, upper;
8778 int anglesNormal = (
simParams->extraBondsCosAngles ? 0 : 1 );
8779 #ifndef MEM_OPT_VERSION 8792 NAMD_die(
"NO EXTRA BONDS FILES SPECIFIED");
8795 for ( ; file; file = file->
next ) {
8796 FILE *f = fopen(file->
data,
"r");
8798 sprintf(err_msg,
"UNABLE TO OPEN EXTRA BONDS FILE %s", file->
data);
8810 if (ret_code!=0)
break;
8813 sscanf(buffer,
"%s",type);
8815 #define CHECKATOMID(ATOMID) if ( ATOMID < 0 || ATOMID >= numAtoms ) badatom = 1; 8819 if ( ! strncasecmp(type,
"bond",4) ) {
8820 if ( sscanf(buffer,
"%s %d %d %f %f %s",
8821 type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1;
8827 #ifndef MEM_OPT_VERSION 8835 tmpv.
k = k; tmpv.
x0 = ref;
8837 }
else if ( ! strncasecmp(type,
"wall",4) ) {
8840 if ( sscanf(buffer,
"%s %d %d %f %f %f %s",
8841 type, &a1, &a2, &k, &ref, &upper, err_msg) != 6 ) badline = 1;
8842 else if (upper < ref) badline = 1;
8848 #ifndef MEM_OPT_VERSION 8856 tmpv.
k = k; tmpv.
x0 = ref; tmpv.
x1 = upper;
8858 }
else if ( ! strncasecmp(type,
"angle",4) ) {
8859 if ( sscanf(buffer,
"%s %d %d %d %f %f %s",
8860 type, &a1, &a2, &a3, &k, &ref, err_msg) != 6 ) badline = 1;
8866 #ifndef MEM_OPT_VERSION 8874 tmpv.
k = k; tmpv.
theta0 = ref / 180. *
PI;
8876 tmpv.
normal = anglesNormal;
8879 }
else if ( ! strncasecmp(type,
"dihedral",4) ) {
8881 int ret = 1 + sscanf(buffer,
"%s %d %d %d %d %f %f %s",
8882 type, &a1, &a2, &a3, &a4, &k, &ref, err_msg);
8884 ret = sscanf(buffer,
"%s %d %d %d %d %f %d %f %s",
8885 type, &a1, &a2, &a3, &a4, &k, &n, &ref, err_msg);
8887 if ( ret != 8 ) badline = 1;
8894 #ifndef MEM_OPT_VERSION 8905 }
else if ( ! strncasecmp(type,
"improper",4) ) {
8907 int ret = 1 + sscanf(buffer,
"%s %d %d %d %d %f %f %s",
8908 type, &a1, &a2, &a3, &a4, &k, &ref, err_msg);
8910 ret = sscanf(buffer,
"%s %d %d %d %d %f %d %f %s",
8911 type, &a1, &a2, &a3, &a4, &k, &n, &ref, err_msg);
8913 if ( ret != 8 ) badline = 1;
8920 #ifndef MEM_OPT_VERSION 8931 }
else if ( ! strncasecmp(type,
"#",1) ) {
8938 sprintf(err_msg,
"BAD LINE IN EXTRA BONDS FILE %s: %s",
8939 file->
data, buffer);
8943 sprintf(err_msg,
"BAD ATOM ID IN EXTRA BONDS FILE %s: %s",
8944 file->
data, buffer);
8953 if ( extraNumBonds ) {
8954 iout <<
iINFO <<
"READ " << extraNumBonds <<
" EXTRA BONDS\n" <<
endi;
8956 #ifndef MEM_OPT_VERSION 8957 Bond *newbonds =
new Bond[numBonds+extraNumBonds];
8958 memcpy(newbonds, this->bonds, numBonds*
sizeof(
Bond));
8959 memcpy(newbonds+numBonds, bonds.
begin(), extraNumBonds*
sizeof(
Bond));
8960 delete [] this->bonds;
8961 this->bonds = newbonds;
8962 numBonds += extraNumBonds;
8977 if ( extraNumAngles ) {
8978 iout <<
iINFO <<
"READ " << extraNumAngles <<
" EXTRA ANGLES\n" <<
endi;
8979 if ( anglesNormal ) {
8980 iout <<
iINFO <<
"EXTRA ANGLES ARE NORMAL HARMONIC\n" <<
endi;
8981 }
else if (
simParams->extraBondsCosAnglesSetByUser ) {
8984 iout <<
iWARN <<
"EXTRA ANGLES ARE COSINE-BASED BY DEFAULT TO MATCH PREVIOUS VERSIONS\n";
8985 iout <<
iWARN <<
"FOR NORMAL HARMONIC EXTRA ANGLES SET extraBondsCosAngles off\n" <<
endi;
8987 #ifndef MEM_OPT_VERSION 8988 Angle *newangles =
new Angle[numAngles+extraNumAngles];
8989 memcpy(newangles, this->angles, numAngles*
sizeof(
Angle));
8990 memcpy(newangles+numAngles, angles.
begin(), extraNumAngles*
sizeof(
Angle));
8991 delete [] this->angles;
8992 this->angles = newangles;
8993 numAngles += extraNumAngles;
9008 if ( extraNumDihedrals ) {
9009 iout <<
iINFO <<
"READ " << extraNumDihedrals <<
" EXTRA DIHEDRALS\n" <<
endi;
9010 #ifndef MEM_OPT_VERSION 9012 memcpy(newdihedrals, this->dihedrals, numDihedrals*
sizeof(
Dihedral));
9013 memcpy(newdihedrals+numDihedrals, dihedrals.
begin(), extraNumDihedrals*
sizeof(
Dihedral));
9014 delete [] this->dihedrals;
9015 this->dihedrals = newdihedrals;
9016 numDihedrals += extraNumDihedrals;
9031 if ( extraNumImpropers ) {
9032 iout <<
iINFO <<
"READ " << extraNumImpropers <<
" EXTRA IMPROPERS\n" <<
endi;
9033 #ifndef MEM_OPT_VERSION 9035 memcpy(newimpropers, this->impropers, numImpropers*
sizeof(
Improper));
9036 memcpy(newimpropers+numImpropers, impropers.
begin(), extraNumImpropers*
sizeof(
Improper));
9037 delete [] this->impropers;
9038 this->impropers = newimpropers;
9039 numImpropers += extraNumImpropers;
9073 PDB *initial_pdb,
char *cwd,
9074 const char *simmethod) {
9084 if (alchfile == NULL) {
9085 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, alchfile required.");
9087 strcpy(filename,
"coordinate pdb file (default)");
9090 if (alchfile->
next != NULL) {
9091 char *new_err_msg =
new char[24 + strlen(simmethod) + 26];
9092 sprintf(new_err_msg,
"Multiple definitions of %sFile in configuration file",simmethod);
9096 if ((cwd == NULL) || (alchfile->
data[0] ==
'/')) {
9097 strcpy(filename, alchfile->
data);
9100 strcpy(filename, cwd);
9101 strcat(filename, alchfile->
data);
9104 bPDB =
new PDB(filename);
9106 NAMD_die(
"Memory allocation failed in Molecule:build_fep_flags");
9110 char *new_err_msg =
new char[19 + strlen(simmethod) + 38];
9111 sprintf(new_err_msg,
"Number of atoms in %sFile PDB does not match coordinate PDB",simmethod);
9119 if (alchcol == NULL) {
9123 if (alchcol->
next != NULL) {
9124 char *new_err_msg =
new char[24 + strlen(simmethod) + 35];
9125 sprintf(new_err_msg,
"Multiple definitions of %s parameter column in config file",simmethod);
9129 if (strcasecmp(alchcol->
data,
"X") == 0) {
9132 else if (strcasecmp(alchcol->
data,
"Y") == 0) {
9135 else if (strcasecmp(alchcol->
data,
"Z") == 0) {
9138 else if (strcasecmp(alchcol->
data,
"O") == 0) {
9141 else if (strcasecmp(alchcol->
data,
"B") == 0) {
9145 NAMD_die(
"alchcol must have value of X, Y, Z, O or B");
9149 iout <<
iINFO <<
"To read " << simmethod <<
"data from file: " << filename
9151 iout <<
iINFO <<
"To read " << simmethod <<
"flag data from column: " << bcol
9155 fepAtomFlags =
new unsigned char[numAtoms];
9157 if (fepAtomFlags == NULL) {
9158 NAMD_die(
"Memory allocation failed in Molecule::build_fep_params()");
9161 double lesMassFactor = 1.0;
9163 lesMassFactor = 1.0 /
simParams->lesFactor;
9167 for (i = 0; i < numAtoms; i++) {
9171 bval = (bPDB->
atom(i))->xcoor();
9174 bval = (bPDB->
atom(i))->ycoor();
9177 bval = (bPDB->
atom(i))->zcoor();
9180 bval = (bPDB->
atom(i))->occupancy();
9183 bval = (bPDB->
atom(i))->temperaturefactor();
9189 if ( bval == (
int) bval && bval > 0 ) {
9191 NAMD_die(
"LES flag must be less than or equal to lesFactor.");
9192 fepAtomFlags[i] = (int) bval;
9193 #ifdef MEM_OPT_VERSION 9194 Real newMass = atomMassPool[eachAtomMass[i]]*lesMassFactor;
9195 eachAtomMass[i] = insert_new_mass(newMass);
9197 atoms[i].mass *= lesMassFactor;
9202 fepAtomFlags[i] = 0;
9206 fepAtomFlags[i] = 3;
9208 }
else if (bval == 1.0) {
9211 }
else if (bval == -1.0) {
9212 fepAtomFlags[i] = 2;
9214 }
else if (bval == -2.0) {
9215 fepAtomFlags[i] = 4;
9218 fepAtomFlags[i] = 0;
9220 }
else if (
simParams->pairInteractionOn) {
9221 if (bval ==
simParams->pairInteractionGroup1) {
9222 fepAtomFlags[i] = 1;
9224 }
else if (bval ==
simParams->pairInteractionGroup2) {
9225 fepAtomFlags[i] = 2;
9228 fepAtomFlags[i] = 0;
9230 }
else if (
simParams->pressureProfileAtomTypes > 1) {
9231 fepAtomFlags[i] = (int) bval;
9233 #ifdef OPENATOM_VERSION 9237 fepAtomFlags[i] = bval;
9240 fepAtomFlags[i] = 0;
9243 #endif //OPENATOM_VERSION 9247 if (alchfile != NULL) {
9271 if (ssfile == NULL) {
9272 if ( ! initial_pdb ) {
9273 NAMD_die(
"Initial PDB file unavailable, soluteScalingFile required.");
9276 strcpy(filename,
"coordinate PDB file (default)");
9279 if (ssfile->
next != NULL) {
9280 NAMD_die(
"Multiple definitions of soluteScalingFile in configuration file");
9283 if ((cwd == NULL) || (ssfile->
data[0] ==
'/')) {
9284 strcpy(filename, ssfile->
data);
9287 strcpy(filename, cwd);
9288 strcat(filename, ssfile->
data);
9291 bPDB =
new PDB(filename);
9293 NAMD_die(
"Memory allocation failed in Molecule::build_ss_flags");
9297 NAMD_die(
"Number of atoms in soluteScalingFile PDB does not match coordinate PDB");
9301 if (sscol == NULL) {
9305 if (sscol->
next != NULL) {
9306 NAMD_die(
"Multiple definitions of soluteScalingCol value in config file");
9309 if (strcasecmp(sscol->
data,
"X") == 0) {
9312 else if (strcasecmp(sscol->
data,
"Y") == 0) {
9315 else if (strcasecmp(sscol->
data,
"Z") == 0) {
9318 else if (strcasecmp(sscol->
data,
"O") == 0) {
9321 else if (strcasecmp(sscol->
data,
"B") == 0) {
9325 NAMD_die(
"soluteScalingCol must have value of X, Y, Z, O or B");
9329 iout <<
iINFO <<
"Reading solute scaling data from file: " 9330 << filename <<
"\n" <<
endi;
9331 iout <<
iINFO <<
"Reading solute scaling flags from column: " 9332 << bcol <<
"\n" <<
endi;
9334 ssAtomFlags =
new unsigned char[numAtoms];
9335 ss_index =
new int[numAtoms];
9337 if (ssAtomFlags == NULL || ss_index == NULL) {
9338 NAMD_die(
"Memory allocation failed in Molecule::build_ss_params()");
9342 for (i = 0; i < numAtoms; i++) {
9345 bval = (bPDB->
atom(i))->xcoor();
9348 bval = (bPDB->
atom(i))->ycoor();
9351 bval = (bPDB->
atom(i))->zcoor();
9354 bval = (bPDB->
atom(i))->occupancy();
9357 bval = (bPDB->
atom(i))->temperaturefactor();
9363 ss_index[num_ss] = i;
9372 if (ssfile != NULL) {
9381 int *numAtomsByLjType =
new int[LJtypecount];
9385 ss_vdw_type =
new int[LJtypecount];
9388 for (i = 0; i < LJtypecount; i++) {
9389 numAtomsByLjType[i] = 0;
9394 for (i = 0; i < num_ss; i++) {
9395 numAtomsByLjType[atoms[ss_index[i]].vdw_type]++;
9399 ss_num_vdw_params = 0;
9400 for (i = 0; i < LJtypecount; i++) {
9402 if (numAtomsByLjType[i] != 0) {
9405 ss_vdw_type[ss_num_vdw_params] = i;
9408 ss_num_vdw_params++;
9412 for (i = 0; i < num_ss; i++) {
9414 for (j = 0; j < ss_num_vdw_params; j++) {
9416 if (atoms[ss_index[i]].vdw_type == ss_vdw_type[j]) {
9419 atoms[ss_index[i]].vdw_type = LJtypecount + j;
9424 delete [] numAtomsByLjType;
9438 #ifndef MEM_OPT_VERSION 9442 FILE *alch_unpert_bond_file;
9445 if ((alch_unpert_bond_file =
Fopen(alch_fname,
"r")) == NULL) {
9446 sprintf(err_msg,
"UNABLE TO OPEN ALCH UNPERTBURBED BOND FILE %s", alch_fname);
9456 NAMD_die(
"DID NOT FIND NBOND AFTER ATOM LIST IN ALCH UNPERT PSF");
9460 sscanf(buffer,
"%d", &num_alch_unpert_Bonds);
9462 read_alch_unpert_bonds(alch_unpert_bond_file);
9471 NAMD_die(
"DID NOT FIND NTHETA AFTER BOND LIST IN ALCH UNPERT PSF");
9475 sscanf(buffer,
"%d", &num_alch_unpert_Angles);
9477 read_alch_unpert_angles(alch_unpert_bond_file);
9488 NAMD_die(
"DID NOT FIND NPHI AFTER ANGLE LIST IN ALCH UNPERT PSF");
9492 sscanf(buffer,
"%d", &num_alch_unpert_Dihedrals);
9494 read_alch_unpert_dihedrals(alch_unpert_bond_file);
9496 Fclose(alch_unpert_bond_file);
9502 suspiciousAlchBonds = 0;
9503 for (
int i = 0; i < numBonds; i++) {
9504 int part1 = fepAtomFlags[bonds[i].atom1];
9505 int part2 = fepAtomFlags[bonds[i].atom2];
9506 if ((part1 == 1 || part2 == 1 ) &&
9507 (part1 == 2 || part2 == 2 )) {
9509 suspiciousAlchBonds++;
9514 Angle *nonalchAngles;
9515 nonalchAngles =
new Angle[numAngles];
9516 int nonalchAngleCount = 0;
9517 alchDroppedAngles = 0;
9518 for (
int i = 0; i < numAngles; i++) {
9519 int part1 = fepAtomFlags[angles[i].atom1];
9520 int part2 = fepAtomFlags[angles[i].atom2];
9521 int part3 = fepAtomFlags[angles[i].atom3];
9522 if ((part1 == 1 || part2 == 1 || part3 == 1) &&
9523 (part1 == 2 || part2 == 2 || part3 == 2)) {
9525 alchDroppedAngles++;
9528 if ( angles[i].angle_type == -1 ) {
9531 "MISSING PARAMETERS FOR ANGLE %i %i %i PARTITIONS %i %i %i\n",
9532 angles[i].atom1+1, angles[i].atom2+1, angles[i].atom3+1,
9533 part1, part2, part3);
9536 nonalchAngles[nonalchAngleCount++] = angles[i];
9539 numAngles = nonalchAngleCount;
9541 angles =
new Angle[numAngles];
9542 for (
int i = 0; i < nonalchAngleCount; i++) {
9543 angles[i]=nonalchAngles[i];
9545 delete [] nonalchAngles;
9550 nonalchDihedrals =
new Dihedral[numDihedrals];
9551 int nonalchDihedralCount = 0;
9552 alchDroppedDihedrals = 0;
9553 for (
int i = 0; i < numDihedrals; i++) {
9554 int part1 = fepAtomFlags[dihedrals[i].atom1];
9555 int part2 = fepAtomFlags[dihedrals[i].atom2];
9556 int part3 = fepAtomFlags[dihedrals[i].atom3];
9557 int part4 = fepAtomFlags[dihedrals[i].atom4];
9558 if ((part1 == 1 || part2 == 1 || part3 == 1 || part4 == 1) &&
9559 (part1 == 2 || part2 == 2 || part3 == 2 || part4 == 2)) {
9561 alchDroppedDihedrals++;
9564 if ( dihedrals[i].dihedral_type == -1 ) {
9567 "MISSING PARAMETERS FOR DIHEDRAL %i %i %i %i PARTITIONS %i %i %i %i\n",
9568 dihedrals[i].atom1+1, dihedrals[i].atom2+1,
9569 dihedrals[i].atom3+1, dihedrals[i].atom4+1,
9570 part1, part2, part3, part4);
9573 nonalchDihedrals[nonalchDihedralCount++] = dihedrals[i];
9576 numDihedrals = nonalchDihedralCount;
9577 delete [] dihedrals;
9578 dihedrals =
new Dihedral[numDihedrals];
9579 for (
int i = 0; i < numDihedrals; i++) {
9580 dihedrals[i]=nonalchDihedrals[i];
9582 delete [] nonalchDihedrals;
9586 nonalchImpropers =
new Improper[numImpropers];
9587 int nonalchImproperCount = 0;
9588 alchDroppedImpropers = 0;
9589 for (
int i = 0; i < numImpropers; i++) {
9590 int part1 = fepAtomFlags[impropers[i].atom1];
9591 int part2 = fepAtomFlags[impropers[i].atom2];
9592 int part3 = fepAtomFlags[impropers[i].atom3];
9593 int part4 = fepAtomFlags[impropers[i].atom4];
9594 if ((part1 == 1 || part2 == 1 || part3 == 1 || part4 == 1) &&
9595 (part1 == 2 || part2 == 2 || part3 == 2 || part4 == 2)) {
9597 alchDroppedImpropers++;
9600 nonalchImpropers[nonalchImproperCount++] = impropers[i];
9603 numImpropers = nonalchImproperCount;
9604 delete [] impropers;
9605 impropers =
new Improper[numImpropers];
9606 for (
int i = 0; i < numImpropers; i++) {
9607 impropers[i]=nonalchImpropers[i];
9609 delete [] nonalchImpropers;
9630 if (fixedfile == NULL) {
9631 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, excludeFromPressureFile required.");
9634 if (fixedfile->
next != NULL) {
9635 NAMD_die(
"Multiple definitions of excluded pressure atoms PDB file in configuration file");
9638 if ( (cwd == NULL) || (fixedfile->
data[0] ==
'/') ) {
9639 strcpy(filename, fixedfile->
data);
9641 strcpy(filename, cwd);
9642 strcat(filename, fixedfile->
data);
9644 bPDB =
new PDB(filename);
9645 if ( bPDB == NULL ) {
9646 NAMD_die(
"Memory allocation failed in Molecule::build_exPressure_atoms");
9650 NAMD_die(
"Number of atoms in excludedPressure atoms PDB doesn't match coordinate PDB");
9659 if (fixedcol == NULL) {
9662 if (fixedcol->
next != NULL) {
9663 NAMD_die(
"Multiple definitions of excludedPressure atoms column in config file");
9666 if (strcasecmp(fixedcol->
data,
"X") == 0) {
9668 }
else if (strcasecmp(fixedcol->
data,
"Y") == 0) {
9670 }
else if (strcasecmp(fixedcol->
data,
"Z") == 0) {
9672 }
else if (strcasecmp(fixedcol->
data,
"O") == 0) {
9674 }
else if (strcasecmp(fixedcol->
data,
"B") == 0) {
9677 NAMD_die(
"excludedPressureFileCol must have value of X, Y, Z, O, or B");
9682 exPressureAtomFlags =
new int32[numAtoms];
9684 if (exPressureAtomFlags == NULL) {
9685 NAMD_die(
"memory allocation failed in Molecule::build_fixed_atoms()");
9688 numExPressureAtoms = 0;
9691 for (i=0; i<numAtoms; i++) {
9694 case 1: bval = (bPDB->
atom(i))->xcoor();
break;
9695 case 2: bval = (bPDB->
atom(i))->ycoor();
break;
9696 case 3: bval = (bPDB->
atom(i))->zcoor();
break;
9697 case 4: bval = (bPDB->
atom(i))->occupancy();
break;
9698 case 5: bval = (bPDB->
atom(i))->temperaturefactor();
break;
9703 exPressureAtomFlags[i] = 1;
9704 numExPressureAtoms++;
9706 exPressureAtomFlags[i] = 0;
9709 if (fixedfile != NULL)
9712 iout <<
iINFO <<
"Got " << numExPressureAtoms <<
" excluded pressure atoms." 9722 return ((atoms[anum].status &
DrudeAtom) != 0);
9732 return ((atoms[anum].status &
OxygenAtom) != 0);
9737 return (hydrogenGroup[atoms[anum].hydrogenList].isGP);
9742 return (hydrogenGroup[atoms[anum].hydrogenList].waterVal == 2);
9747 return (hydrogenGroup[atoms[anum].hydrogenList].atomsInGroup);
9754 return atoms[anum].partner;
9758 if ( n != numAtoms )
9759 NAMD_die(
"Incorrect number of atoms in Molecule::reloadCharges().");
9761 #ifdef MEM_OPT_VERSION 9762 delete [] atomChargePool;
9763 vector<Real> tmpCharges;
9764 for(
int i=0; i<numAtoms; i++){
9768 for(
int j=0; j<tmpCharges.size();j++){
9769 if(tmpCharges[j] == charge[i]){
9775 tmpCharges.push_back(charge[i]);
9776 foundIdx = tmpCharges.size()-1;
9778 eachAtomCharge[i] = (
Index)foundIdx;
9780 chargePoolSize = tmpCharges.size();
9781 atomChargePool =
new Real[chargePoolSize];
9782 for(
int i=0; i<chargePoolSize; i++)
9783 atomChargePool[i] = tmpCharges[i];
9785 for(
int i=0; i<n; ++i ) atoms[i].charge = charge[i];
9797 #ifndef MEM_OPT_VERSION 9798 if(std::filesystem::path(dcdSelectionInputFile).extension() ==
".idx")
9801 IndexFile dcdSelectionInputIdx(dcdSelectionInputFile);
9802 std::vector<uint32> indexVec = dcdSelectionInputIdx.getAllElements();
9803 dcdSelectionParams[index].size = dcdSelectionInputIdx.size();
9804 int bitmask = 1 << index;
9805 for(
int i=0; i<indexVec.size(); ++i)
9806 { atoms[indexVec[i]].flags.dcdSelection |= bitmask; }
9810 PDB dcdSelectionInputFilePdb(dcdSelectionInputFile);
9812 int numDcdSelectionAtoms = dcdSelectionInputFilePdb.num_atoms();
9813 dcdSelectionParams[index].size = numDcdSelectionAtoms;
9814 if (numDcdSelectionAtoms < 1)
9815 NAMD_die(
"No atoms found in dcdSelectionInputFile\n");
9820 if (numDcdSelectionAtoms != numAtoms)
9821 NAMD_die(
"The number of atoms in dcdSelectionInputFile must be equal to the total number of atoms in the structure!");
9822 int bitmask = 1 << index;
9823 for (
int i=0; i<numAtoms; i++) {
9824 PDBAtom *atom = dcdSelectionInputFilePdb.atom(i);
9827 atoms[i].flags.dcdSelection |= bitmask;
9832 if(std::filesystem::path(dcdSelectionInputFile).extension() ==
".idx")
9835 IndexFile dcdSelectionInputIdx(dcdSelectionInputFile);
9836 dcdSelectionParams[index].size=dcdSelectionInputIdx.size();
9840 dcdSelectionParams[index].tag = index+1;
9847 #ifndef MEM_OPT_VERSION 9848 dcdSelectionParams[index].dcdSelectionIndexReverse.resize(numAtoms);
9850 int bitmask = 1 << index;
9852 for (
int i=0; i<numAtoms; i++) {
9853 if(atoms[i].flags.dcdSelection & bitmask)
9855 dcdSelectionParams[index].dcdSelectionIndex.push_back(i);
9856 dcdSelectionParams[index].dcdSelectionIndexReverse[i] = count++;
9860 dcdSelectionParams[index].dcdSelectionIndexReverse[i] = -1;
9863 dcdSelectionParams[index].size = count;
9871 dcdSelectionParams[index].frequency = freq;
9876 if(
auto keyIt=dcdSelectionKeyMap.find(std::string(keystr)); keyIt!=dcdSelectionKeyMap.end())
9878 return(keyIt->second);
9886 uint16 key=find_dcd_selection_index(keystr);
9893 key = dcdSelectionKeyMap.size();
9897 sprintf(errmsg,
"Key %s beyond Max (16) DCD Selections\n", keystr);
9900 dcdSelectionKeyMap[std::string(keystr)] = key;
9916 char *keystr =
new char[256];
9917 char *valstr =
new char[256];
9920 int dcdSelectionCounter[3]={0,0,0};
9921 StringList *dcdSelectionInputFileList = config->
find(
"dcdselectioninputfile");
9922 for(;dcdSelectionInputFileList; dcdSelectionInputFileList = dcdSelectionInputFileList->
next)
9924 sscanf(dcdSelectionInputFileList->
data,
"%80s %255s",keystr,valstr);
9925 key = find_or_create_dcd_selection_index(keystr);
9926 dcdSelectionCounter[0]++;
9927 build_dcd_selection_list_pdb(key, valstr);
9928 iout <<
iINFO <<
"Dcdselection tag "<< keystr <<
" Input file " << valstr <<
" index "<< key <<
"\n";
9930 StringList *dcdSelectionFileList = config->
find(
"dcdselectionfile");
9931 for(;dcdSelectionFileList; dcdSelectionFileList = dcdSelectionFileList->
next)
9934 sscanf(dcdSelectionFileList->
data,
"%80s %255s",keystr,valstr);
9935 key = find_dcd_selection_index(keystr);
9938 add_dcd_selection_file(key, valstr);
9939 dcdSelectionCounter[1]++;
9943 sprintf(errmsg,
"Keystring %s has no dcd input file match\n", keystr);
9946 iout <<
iINFO <<
"Dcdselection tag "<< keystr <<
" Output file " << valstr<<
" index "<< key <<
"\n";
9948 StringList *dcdSelectionFreqList = config->
find(
"dcdselectionfreq");
9949 for(;dcdSelectionFreqList; dcdSelectionFreqList = dcdSelectionFreqList->
next)
9952 sscanf(dcdSelectionFreqList->
data,
"%s %d",keystr,&freq);
9953 key = find_dcd_selection_index(keystr);
9956 add_dcd_selection_freq(key, freq);
9957 dcdSelectionCounter[2]++;
9961 sprintf(errmsg,
"Keystring %s has no dcd input file match\n", keystr);
9964 iout <<
iINFO <<
"Dcd Selection tag "<< keystr <<
" freq " << freq<<
" index "<< key <<
"\n";
9968 if((dcdSelectionCounter[0]!=dcdSelectionCounter[1] && CmiNumPartitions()==1) || dcdSelectionCounter[0]!=dcdSelectionCounter[2] || dcdSelectionCounter[0]>15)
9971 sprintf(errmsg,
"Invalid DCD Selection configuration\n");
9977 #ifndef MEM_OPT_VERSION 9984 void Molecule::build_atom_status(
void) {
9987 int numDrudeWaters = 0;
9990 numLonepairs = numZeroMassAtoms;
9992 iout <<
iWARN <<
"CORRECTION OF ZERO MASS ATOMS TURNED OFF " 9993 "BECAUSE LONE PAIRS ARE USED\n" <<
endi;
9997 if (is_lonepairs_psf && numLonepairs != numLphosts) {
9998 NAMD_die(
"must have same number of LP hosts as lone pairs");
10000 }
else if (numZeroMassAtoms) {
10001 for (i=0; i < numAtoms; i++) {
10002 if ( atoms[i].mass < 0.001 ) atoms[i].mass = 0.001;
10004 if ( ! CkMyPe() ) {
10005 iout <<
iWARN <<
"FOUND " << numZeroMassAtoms <<
10006 " ATOMS WITH ZERO OR NEGATIVE MASSES! CHANGED TO 0.001\n" <<
endi;
10011 hydrogenGroup.resize(numAtoms);
10013 for (i=0; i < numAtoms; i++) {
10014 atoms[i].partner = (-1);
10024 int hhbondcount = 0;
10025 for (i=0; i < numRealBonds; i++) {
10026 a1 = bonds[i].atom1;
10027 a2 = bonds[i].atom2;
10028 if (is_hydrogen(a1) && is_hydrogen(a2)) {
10031 atoms[a1].partner = a2;
10032 atoms[a2].partner = a1;
10040 if ( hhbondcount && ! CkMyPe() ) {
10041 iout <<
iWARN <<
"Found " << hhbondcount <<
" H-H bonds.\n" <<
endi;
10046 for (i=0; i < numRealBonds; i++) {
10047 a1 = bonds[i].atom1;
10048 a2 = bonds[i].atom2;
10049 if (is_hydrogen(a1)) {
10050 if (is_hydrogen(a2))
continue;
10051 atoms[a1].partner = a2;
10057 if (is_oxygen(a2)) hg[a2].waterVal++;
10059 if (is_hydrogen(a2)) {
10060 atoms[a2].partner = a1;
10066 if (is_oxygen(a1)) hg[a1].waterVal++;
10072 atoms[a1].partner = a2;
10079 atoms[a2].partner = a1;
10087 else if ( is_lonepairs_psf || is_drude_psf ) {
10088 if (is_lp(a1) || is_drude(a1)) {
10089 if (is_hydrogen(a2) || is_lp(a2) || is_drude(a2)) {
10091 sprintf(msg,
"%s particle %d is bonded to non-parent atom %d",
10092 (is_lp(a1) ?
"Lone pair" :
"Drude"), a1+1, a2+1);
10095 atoms[a1].partner = a2;
10101 else if (is_lp(a2) || is_drude(a2)) {
10102 if (is_hydrogen(a1) || is_lp(a1) || is_drude(a1)) {
10104 sprintf(msg,
"%s particle %d is bonded to non-parent atom %d",
10105 (is_lp(a2) ?
"Lone pair" :
"Drude"), a2+1, a1+1);
10108 atoms[a2].partner = a1;
10120 for(i=0; i<numAtoms; i++) {
10121 if ( ! hg[hg[i].GPID].isGP ) {
10123 sprintf(msg,
"child atom %d bonded only to child H atoms",i+1);
10126 if ( hg[i].isGP && is_hydrogen(i) ) {
10127 if ( hg[i].GPID == i )
continue;
10129 if ( is_hydrogen(hg[i].GPID) && hg[hg[i].GPID].GPID != i ) {
10131 sprintf(msg,
"H atom %d bonded only to child H atoms",i+1);
10137 if ( hg[i].atomsInGroup != 2 ) {
10139 sprintf(msg,
"H atom %d bonded to multiple H atoms",i+1);
10144 if ( hGPcount && ! CkMyPe() ) {
10145 iout <<
iWARN <<
"Found " << hGPcount <<
" H-H molecules.\n" <<
endi;
10149 for (i=0; i<numAtoms; ++i) {
10150 if ( hg[i].isGP ) hg[i].
GPID = i;
10158 for (i=0; i<numLphosts; ++i) {
10159 int a1 = lphosts[i].atom1;
10160 int a2 = lphosts[i].atom2;
10161 int a3 = lphosts[i].atom3;
10162 int a4 = lphosts[i].atom4;
10163 int m1 = hg[a1].
MPID;
10164 while ( hg[m1].MPID != m1 ) m1 = hg[m1].
MPID;
10165 int m2 = hg[a2].
MPID;
10166 while ( hg[m2].MPID != m2 ) m2 = hg[m2].
MPID;
10167 int m3 = hg[a3].
MPID;
10168 while ( hg[m3].MPID != m3 ) m3 = hg[m3].
MPID;
10169 int m4 = hg[a4].
MPID;
10170 while ( hg[m4].MPID != m4 ) m4 = hg[m4].
MPID;
10172 if ( m2 < mp ) mp = m2;
10173 if ( m3 < mp ) mp = m3;
10174 if ( m4 < mp ) mp = m4;
10182 for (i=0; i<numAtoms; ++i) {
10183 int mp = hg[i].
MPID;
10184 if ( hg[mp].MPID != mp ) {
10189 if ( allok )
break;
10191 for (i=0; i<numAtoms; ++i) {
10195 for (i=0; i<numAtoms; ++i) {
10201 for (i=0; i<numAtoms; i++) {
10212 numHydrogenGroups = 0;
10213 maxHydrogenGroupSize = 0;
10214 numMigrationGroups = 0;
10215 maxMigrationGroupSize = 0;
10216 for(i=0; i<numAtoms; i++)
10219 ++numMigrationGroups;
10221 if ( mgs > maxMigrationGroupSize ) maxMigrationGroupSize = mgs;
10224 ++numHydrogenGroups;
10226 if ( hgs > maxHydrogenGroupSize ) maxHydrogenGroupSize = hgs;
10230 hydrogenGroup.sort();
10235 for(i=0; i<numAtoms; ++i, --hgs) {
10237 if ( hg[i].isGP ) {
10239 parentid = hg[i].
atomID;
10242 sprintf(buff,
"Atom %d has bad hydrogen group size. " 10243 "Check for duplicate bonds.", parentid+1);
10247 if ( hg[i].isGP ) {
10249 sprintf(buff,
"Atom %d has bad hydrogen group size. " 10250 "Check for duplicate bonds.", parentid+1);
10258 for(i=0; i<numAtoms; ++i, --mgs) {
10260 if ( hg[i].isMP ) {
10262 parentid = hg[i].
atomID;
10265 sprintf(buff,
"Atom %d has bad migration group size.", parentid+1);
10269 if ( hg[i].isMP ) {
10271 sprintf(buff,
"Atom %d has bad migration group size.", parentid+1);
10279 for(i=0; i<numAtoms; i++) {
10280 atoms[hydrogenGroup[i].atomID].hydrogenList = i;
10286 for (i = 0; i < numAtoms; i++) {
10287 if (is_water(hg[i].atomID) && hg[i].isGP) {
10289 || ! is_drude(hg[i+1].atomID)
10290 || ! is_lp(hg[i+2].atomID)
10291 || ! is_hydrogen(hg[i+3].atomID)
10292 || ! is_hydrogen(hg[i+4].atomID) ) {
10294 sprintf(msg,
"Drude water molecule from HydrogenGroup i=%d " 10295 "starting at atom %d is not sorted\n", i, hg[i].atomID+1);
10302 else if (is_drude(hg[i].atomID)) {
10303 if (i < 1 || hg[i-1].atomID != hg[i].GPID) {
10305 sprintf(msg,
"Drude particle from HydrogenGroup i=%d must " 10306 "immediately follow its parent atom %d\n", i, hg[i].GPID+1);
10311 else if (is_lp(hg[i].atomID)) {
10313 sprintf(msg,
"Drude lonepair from HydrogenGroup i=%d " 10314 "at particle %d is NOT from water - unsupported\n",
10315 i, hg[i].atomID+1);
10325 for(i=0; i<numAtoms; i++)
10326 iout << i <<
" atomID=" << hydrogenGroup[i].atomID
10327 <<
" isGP=" << hydrogenGroup[i].isGP
10328 <<
" parent=" << hydrogenGroup[i].GPID
10329 <<
" #" << hydrogenGroup[i].atomsInGroup
10330 <<
" waterVal=" << hydrogenGroup[i].waterVal
10331 <<
" partner=" << atoms[i].partner
10332 <<
" hydrogenList=" << atoms[i].hydrogenList
10343 delete [] rigidBondLengths;
10344 rigidBondLengths =
new Real[numAtoms];
10345 if ( ! rigidBondLengths ) {
10346 NAMD_die(
"Memory allocation failed in Molecule::build_atom_status()\n");
10348 for (i=0; i<numAtoms; ++i) rigidBondLengths[i] = 0;
10353 for (i=0; i < numRealBonds; i++) {
10354 a1 = bonds[i].atom1;
10355 a2 = bonds[i].atom2;
10358 if (is_hydrogen(a2)) {
int tmp = a1; a1 = a2; a2 = tmp; }
10359 if (is_hydrogen(a1)) {
10360 if ( is_hydrogen(a2) ) {
10361 if ( ! is_water(a2) ) {
10362 rigidBondLengths[a1] = ( mode ==
RIGID_ALL ? x0 : 0. );
10363 rigidBondLengths[a2] = ( mode ==
RIGID_ALL ? x0 : 0. );
10365 }
else if ( is_water(a2) || mode ==
RIGID_ALL ) {
10366 rigidBondLengths[a1] = x0;
10367 if (is_water(a2)) r_oh = rigidBondLengths[a1];
10369 rigidBondLengths[a1] = 0.;
10374 if (is_lp(a2)) {
int tmp = a1; a1 = a2; a2 = tmp; }
10376 if (! is_water(a2) ) {
10380 sprintf(err_msg,
"ILLEGAL LONE PAIR AT INDEX %i\n" 10381 "LONE PAIRS ARE CURRENTLY ALLOWED ONLY ON WATER MOLECULES\n",
10385 rigidBondLengths[a1] = x0;
10394 int tmp = a1; a1 = a2; a2 = tmp;
10397 if (is_water(a2)) {
10404 sprintf(msg,
"ILLEGAL LONE PAIR AT INDEX %d\n" 10405 "LONE PAIRS ARE CURRENTLY ALLOWED ONLY ON WATER MOLECULES\n",
10415 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
10416 for( ; h_i != h_e; ++h_i ) {
10417 if ( h_i->
isGP ) rigidBondLengths[h_i->
atomID] = 0.;
10421 for (i=0; i < numAngles; i++) {
10422 a2 = angles[i].atom2;
10423 if ( ! is_water(a2) )
continue;
10424 if ( ! is_oxygen(a2) )
continue;
10425 a1 = angles[i].atom1;
10426 if ( ! is_hydrogen(a1) )
continue;
10427 a3 = angles[i].atom3;
10428 if ( ! is_hydrogen(a3) )
continue;
10429 if (is_lp(a2) || is_lp(a1) || is_lp(a3) ||
10430 is_drude(a2) || is_drude(a1) || is_drude(a3))
continue;
10431 if ( rigidBondLengths[a1] != rigidBondLengths[a3] ) {
10432 if (rigidBondLengths[a1] >0.3 && rigidBondLengths[a3] >0.3) {
10433 printf(
"length1: %f length2: %f\n", rigidBondLengths[a1], rigidBondLengths[a3]);
10435 NAMD_die(
"Asymmetric water molecule found??? This can't be right.\n");
10440 rigidBondLengths[a2] = 2. * rigidBondLengths[a1] * sin(0.5*t0);
10441 r_hh = rigidBondLengths[a2];
10445 int numBondWaters = 0;
10446 int numFailedWaters = 0;
10448 for (i=0; i < numRealBonds; i++) {
10449 a1 = bonds[i].atom1;
10450 a2 = bonds[i].atom2;
10451 if ( ! is_hydrogen(a1) )
continue;
10452 if ( ! is_hydrogen(a2) )
continue;
10453 int ma1 = get_mother_atom(a1);
10454 int ma2 = get_mother_atom(a2);
10455 if ( ma1 != ma2 )
continue;
10456 if ( ! is_water(ma1) )
continue;
10457 if ( rigidBondLengths[ma1] != 0. )
continue;
10460 rigidBondLengths[ma1] = x0;
10467 if (r_oh < 0.0 || r_hh < 0.0) {
10469 NAMD_die(
"Failed to find water bond lengths\n");
10471 r_ohc = sqrt(r_oh * r_oh - 0.25 * r_hh * r_hh);
10475 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
10476 for( ; h_i != h_e; ++h_i ) {
10478 rigidBondLengths[h_i->
atomID] == 0. ) {
10479 if ( h_i + 1 == h_e || h_i + 2 == h_e ||
10480 h_i[1].isGP || h_i[2].isGP || h_i->
atomsInGroup != 3 ) {
10481 NAMD_die(
"Abnormal water detected.");
10483 if ( CkNumNodes() > 1 ) {
10484 NAMD_die(
"Unable to determine H-H distance for rigid water because structure has neither H-O-H angle nor H-H bond.");
10490 get_atomtype(btmp.
atom1),
10491 get_atomtype(btmp.
atom2),
10497 rigidBondLengths[h_i->
atomID] = x0;
10504 if ( numBondWaters + numFailedWaters ) {
10505 iout <<
iWARN <<
"Missing angles for " <<
10506 ( numBondWaters + numFailedWaters ) <<
" waters.\n" <<
endi;
10508 if ( numBondWaters ) {
10509 iout <<
iWARN <<
"Obtained H-H distance from bond parameters for " <<
10510 numBondWaters <<
" waters.\n" <<
endi;
10511 iout <<
iWARN <<
"This would not be possible in a multi-process run.\n" <<
endi;
10513 if ( numFailedWaters ) {
10514 iout <<
iERROR <<
"Failed to obtain H-H distance from angles or bonds for " <<
10515 numFailedWaters <<
" waters.\n" <<
endi;
10523 for (i=0; i<numAtoms; ++i) rigidBondLengths[i] *= -1;
10525 for (i=0; i<numAtoms; ++i) {
10526 if ( ! is_water(i) ) rigidBondLengths[i] *= -1;
10532 for (i=0; i<numAtoms; ++i) {
10533 if ( rigidBondLengths[i] > 0. ) ++numRigidBonds;
10554 BigReal LJAvgA, LJAvgB, LJAvgA1, LJAvgB1, LJAvgA2, LJAvgB2;
10555 int numLJsites, numLJsites1, numLJsites2;
10567 const int ss_dim = soluteScalingOn ? ss_num_vdw_params : 0;
10575 if (soluteScalingOn && (ss_vdw_type == NULL)) {
10576 NAMD_bug(
"Solute scaling is used but ss_vdw_type is NULL. " 10577 "It is likely that compute_LJcorrection() is called " 10578 "before build_ss_flags().");
10580 Real A, B, A14, B14;
10581 Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
10582 Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
10583 Real *ATable =
new Real[LJtypecount*LJtypecount];
10584 Real *BTable =
new Real[LJtypecount*LJtypecount];
10585 int useGeom =
simParams->vdwGeometricSigma;
10587 for (
int i = 0; i < LJtypecount; i++) {
10588 for (
int j = 0; j < LJtypecount; j++) {
10589 const int i_type = soluteScalingOn ? ((i >= table_dim_org) ? ss_vdw_type[i-table_dim_org] : i) : i;
10590 const int j_type = soluteScalingOn ? ((j >= table_dim_org) ? ss_vdw_type[j-table_dim_org] : j) : j;
10592 ATable[i*LJtypecount + j] = A;
10593 BTable[i*LJtypecount + j] = B;
10596 params->
get_vdw_params(&sigma_i,&epsilon_i,&sigma_i14,&epsilon_i14,i_type);
10597 params->
get_vdw_params(&sigma_j,&epsilon_j,&sigma_j14,&epsilon_j14,j_type);
10599 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10600 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10601 sigma_ij *= sigma_ij*sigma_ij;
10602 sigma_ij *= sigma_ij;
10603 ATable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10604 BTable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij;
10609 int *numAtomsByLjType =
new int[LJtypecount];
10610 for (
int i=0; i < LJtypecount; i++) {numAtomsByLjType[i]=0;}
10611 for (
int i=0; i < numAtoms; i++) {numAtomsByLjType[atoms[i].vdw_type]++;}
10618 for (
int i=0; i < LJtypecount; i++) {
10619 for (
int j=0; j < LJtypecount; j++) {
10620 A = ATable[i*LJtypecount + j];
10621 B = BTable[i*LJtypecount + j];
10622 if (!A && !B)
continue;
10623 npairs = (numAtomsByLjType[i] - int(i==j))*
BigReal(numAtomsByLjType[j]);
10624 sumOfAs += npairs*A;
10625 sumOfBs += npairs*B;
10627 if (i==j) numLJsites += numAtomsByLjType[i];
10630 delete [] numAtomsByLjType;
10634 LJAvgA = sumOfAs / count;
10635 LJAvgB = sumOfBs / count;
10656 numLJsites1 = numLJsites2 = numLJsites;
10657 int alch_counter = 0;
10658 for (
int i=0; i < numAtoms; ++i) {
10660 if (get_fep_type(i) == 2 || get_fep_type(i) == 4) alchFlagi = -1;
10661 if (get_fep_type(i) == 1 || get_fep_type(i) == 3) alchFlagi = 1;
10663 &A, &B, &A14, &B14)) {
10667 &epsilon_i14, atoms[i].vdw_type);
10669 useGeom ? sqrt(sigma_i*sigma_i) : 0.5*(sigma_i+sigma_i);
10670 BigReal epsilon_ii = sqrt(epsilon_i*epsilon_i);
10672 sigma_ii *= sigma_ii*sigma_ii;
10673 sigma_ii *= sigma_ii;
10674 A = 4.0*sigma_ii*epsilon_ii*sigma_ii;
10675 B = 4.0*sigma_ii*epsilon_ii;
10678 if (alchFlagi == 1) numLJsites2--;
10679 else if (alchFlagi == -1) numLJsites1--;
10681 for (
int j=i+1; j < numAtoms; ++j) {
10683 if (get_fep_type(j) == 2 || get_fep_type(j) == 4) alchFlagj = -1;
10684 if (get_fep_type(j) == 1 || get_fep_type(j) == 3) alchFlagj = 1;
10685 int alchFlagSum = alchFlagi + alchFlagj;
10688 if (alchFlagi == 0 && alchFlagj == 0)
continue;
10691 &A, &B, &A14, &B14)) {
10695 &epsilon_i14, atoms[i].vdw_type);
10697 &epsilon_j14, atoms[j].vdw_type);
10699 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10700 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10702 sigma_ij *= sigma_ij*sigma_ij;
10703 sigma_ij *= sigma_ij;
10704 A = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10705 B = 4.0*sigma_ij*epsilon_ij;
10707 if (!A && !B)
continue;
10712 if ( alchFlagSum > 0 ){
10717 else if ( alchFlagSum < 0 ){
10733 if ( alchFlagi == 1 || alchFlagi == -1 ) alch_counter++;
10734 if ( alch_counter == (numFepInitial + numFepFinal) )
break;
10736 LJAvgA = sumOfAs / count;
10737 LJAvgB = sumOfBs / count;
10739 LJAvgA1 = sumOfAs1 / count1;
10740 LJAvgB1 = sumOfBs1 / count1;
10746 LJAvgA2 = sumOfAs2 / count2;
10747 LJAvgB2 = sumOfBs2 / count2;
10752 if ( ! CkMyPe() ) {
10753 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 10754 <<
"ENERGY AND PRESSURE\n" <<
endi;
10755 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A0 AND B0 COEFFICIENTS " 10756 << LJAvgA2 <<
" AND " << LJAvgB2 <<
"\n" <<
endi;
10757 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A1 AND B1 COEFFICIENTS " 10758 << LJAvgA1 <<
" AND " << LJAvgB1 <<
"\n" <<
endi;
10760 numLJsites = (numLJsites1 + numLJsites2 - numLJsites);
10761 LJAvgA1 *=
BigReal(numLJsites1)*numLJsites1;
10762 LJAvgB1 *=
BigReal(numLJsites1)*numLJsites1;
10763 LJAvgA2 *=
BigReal(numLJsites2)*numLJsites2;
10764 LJAvgB2 *=
BigReal(numLJsites2)*numLJsites2;
10767 LJAvgA1 = LJAvgB1 = LJAvgA2 = LJAvgB2 = 0;
10769 if ( ! CkMyPe() ) {
10770 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 10771 <<
"ENERGY AND PRESSURE\n" <<
endi;
10772 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A AND B COEFFICIENTS " 10773 << LJAvgA <<
" AND " << LJAvgB <<
"\n" <<
endi;
10776 LJAvgA *=
BigReal(numLJsites)*numLJsites;
10777 LJAvgB *=
BigReal(numLJsites)*numLJsites;
10786 BigReal rswitch2 = rswitch*rswitch;
10787 BigReal rswitch3 = rswitch*rswitch2;
10788 BigReal rswitch4 = rswitch2*rswitch2;
10789 BigReal rswitch5 = rswitch2*rswitch3;
10790 BigReal rswitch6 = rswitch3*rswitch3;
10816 BigReal int_U_gofr_A, int_rF_gofr_A, int_U_gofr_B, int_rF_gofr_B;
10835 BigReal rsum3 = (rcut + rswitch)*(rcut + rswitch)*(rcut + rswitch);
10836 int_U_gofr_A = int_rF_gofr_A = (16*
PI*(3*rcut4 + 9*rcut3*rswitch
10837 + 11*rcut2*rswitch2 + 9*rcut*rswitch3
10839 / (315*rcut5*rswitch5*rsum3));
10840 int_U_gofr_B = int_rF_gofr_B = -16*
PI / (3*rsum3);
10867 BigReal lnr = log(rcut/rswitch);
10891 int_U_gofr_A = int_rF_gofr_A =\
10892 16*
PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
10893 int_U_gofr_B = int_rF_gofr_B = -4*
PI*lnr / (rcut3 - rswitch3);
10930 int_rF_gofr_A = 8*
PI / (9*rcut9);
10931 int_rF_gofr_B = -4*
PI / (3*rcut3);
10932 int_U_gofr_A = 2*
PI / (9*rcut9);
10933 int_U_gofr_B = -2*
PI / (3*rcut3);
10936 tail_corr_virial = int_rF_gofr_A*LJAvgA + int_rF_gofr_B*LJAvgB;
10937 tail_corr_ener = int_U_gofr_A*LJAvgA + int_U_gofr_B*LJAvgB;
10939 tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1 -
10941 tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1 -
10943 tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2 -
10945 tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2 -
10964 BigReal LJAvgA, LJAvgB, LJAvgA1, LJAvgB1, LJAvgA2, LJAvgB2;
10965 int numLJsites, numLJsites1, numLJsites2;
10977 const int ss_dim = soluteScalingOn ? ss_num_vdw_params : 0;
10985 if (soluteScalingOn && (ss_vdw_type == NULL)) {
10986 NAMD_bug(
"Solute scaling is used but ss_vdw_type is NULL. " 10987 "It is likely that compute_LJcorrection_alternative() is called " 10988 "before build_ss_flags().");
10990 Real A, B, A14, B14;
10991 Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
10992 Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
10993 Real *ATable =
new Real[LJtypecount*LJtypecount];
10994 Real *BTable =
new Real[LJtypecount*LJtypecount];
10995 int useGeom =
simParams->vdwGeometricSigma;
10997 for (
int i = 0; i < LJtypecount; i++) {
10998 for (
int j = 0; j < LJtypecount; j++) {
10999 const int i_type = soluteScalingOn ? ((i >= table_dim_org) ? ss_vdw_type[i-table_dim_org] : i) : i;
11000 const int j_type = soluteScalingOn ? ((j >= table_dim_org) ? ss_vdw_type[j-table_dim_org] : j) : j;
11002 ATable[i*LJtypecount + j] = A;
11003 BTable[i*LJtypecount + j] = B;
11006 params->
get_vdw_params(&sigma_i,&epsilon_i,&sigma_i14,&epsilon_i14,i_type);
11007 params->
get_vdw_params(&sigma_j,&epsilon_j,&sigma_j14,&epsilon_j14,j_type);
11009 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
11010 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
11011 sigma_ij *= sigma_ij*sigma_ij;
11012 sigma_ij *= sigma_ij;
11013 ATable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij*sigma_ij;
11014 BTable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij;
11019 int *numAtomsByLjType =
new int[LJtypecount];
11020 for (
int i=0; i < LJtypecount; i++) {numAtomsByLjType[i]=0;}
11021 for (
int i=0; i < numAtoms; i++) {numAtomsByLjType[atoms[i].vdw_type]++;}
11026 for (
int i=0; i < LJtypecount; i++) {
11027 for (
int j=0; j < LJtypecount; j++) {
11028 A = ATable[i*LJtypecount + j];
11029 B = BTable[i*LJtypecount + j];
11030 npairs =
BigReal(numAtomsByLjType[i])*
BigReal(numAtomsByLjType[j]);
11031 sumOfAs += npairs*A;
11032 sumOfBs += npairs*B;
11035 delete [] numAtomsByLjType;
11059 int alch_counter = 0;
11060 for (
int i=0; i < numAtoms; ++i) {
11062 if (get_fep_type(i) == 2 || get_fep_type(i) == 4) alchFlagi = -1;
11063 if (get_fep_type(i) == 1 || get_fep_type(i) == 3) alchFlagi = 1;
11065 for (
int j=i; j < numAtoms; ++j) {
11067 if (get_fep_type(j) == 2 || get_fep_type(j) == 4) alchFlagj = -1;
11068 if (get_fep_type(j) == 1 || get_fep_type(j) == 3) alchFlagj = 1;
11069 int alchFlagSum = alchFlagi + alchFlagj;
11072 if (alchFlagi == 0 && alchFlagj == 0)
continue;
11075 &A, &B, &A14, &B14)) {
11079 &epsilon_i14, atoms[i].vdw_type);
11081 &epsilon_j14, atoms[j].vdw_type);
11083 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
11084 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
11086 sigma_ij *= sigma_ij*sigma_ij;
11087 sigma_ij *= sigma_ij;
11088 A = 4.0*sigma_ij*epsilon_ij*sigma_ij;
11089 B = 4.0*sigma_ij*epsilon_ij;
11094 if ( alchFlagSum > 0 ){
11097 }
else if ( alchFlagSum < 0 ){
11109 if ( alchFlagi == 1 || alchFlagi == -1 ) alch_counter++;
11110 if ( alch_counter == (numFepInitial + numFepFinal) )
break;
11115 LJAvgA1 = sumOfAs1;
11116 LJAvgB1 = sumOfBs1;
11117 LJAvgA2 = sumOfAs2;
11118 LJAvgB2 = sumOfBs2;
11120 if ( ! CkMyPe() ) {
11121 BigReal inv_sizeSq = 1.0 / (numAtoms * numAtoms);
11122 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 11123 <<
"ENERGY AND PRESSURE\n" <<
endi;
11124 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A0 AND B0 COEFFICIENTS " 11125 << LJAvgA2*inv_sizeSq <<
" AND " << LJAvgB2*inv_sizeSq <<
"\n" <<
endi;
11126 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A1 AND B1 COEFFICIENTS " 11127 << LJAvgA1*inv_sizeSq <<
" AND " << LJAvgB1*inv_sizeSq <<
"\n" <<
endi;
11130 LJAvgA1 = LJAvgB1 = LJAvgA2 = LJAvgB2 = 0;
11132 if ( ! CkMyPe() ) {
11133 BigReal inv_sizeSq = 1.0 / (numAtoms * numAtoms);
11134 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 11135 <<
"ENERGY AND PRESSURE\n" <<
endi;
11136 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A AND B COEFFICIENTS " 11137 << LJAvgA*inv_sizeSq <<
" AND " << LJAvgB*inv_sizeSq <<
"\n" <<
endi;
11148 BigReal rswitch2 = rswitch*rswitch;
11149 BigReal rswitch3 = rswitch*rswitch2;
11150 BigReal rswitch4 = rswitch2*rswitch2;
11151 BigReal rswitch5 = rswitch2*rswitch3;
11152 BigReal rswitch6 = rswitch3*rswitch3;
11178 BigReal int_U_gofr_A, int_rF_gofr_A, int_U_gofr_B, int_rF_gofr_B;
11197 BigReal rsum3 = (rcut + rswitch)*(rcut + rswitch)*(rcut + rswitch);
11198 int_U_gofr_A = int_rF_gofr_A = (16*
PI*(3*rcut4 + 9*rcut3*rswitch
11199 + 11*rcut2*rswitch2 + 9*rcut*rswitch3
11201 / (315*rcut5*rswitch5*rsum3));
11202 int_U_gofr_B = int_rF_gofr_B = -16*
PI / (3*rsum3);
11229 BigReal lnr = log(rcut/rswitch);
11253 int_U_gofr_A = int_rF_gofr_A =\
11254 16*
PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
11255 int_U_gofr_B = int_rF_gofr_B = -4*
PI*lnr / (rcut3 - rswitch3);
11292 int_rF_gofr_A = 8*
PI / (9*rcut9);
11293 int_rF_gofr_B = -4*
PI / (3*rcut3);
11294 int_U_gofr_A = 2*
PI / (9*rcut9);
11295 int_U_gofr_B = -2*
PI / (3*rcut3);
11298 tail_corr_virial = int_rF_gofr_A*LJAvgA + int_rF_gofr_B*LJAvgB;
11299 tail_corr_ener = int_U_gofr_A*LJAvgA + int_U_gofr_B*LJAvgB;
11301 tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1 -
11303 tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1 -
11305 tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2 -
11307 tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2 -
11319 const BigReal corr = (doti ? 0.0 : tail_corr_ener);
11320 return scl*(corr + vdw_lambda_1*tail_corr_dUdl_1 +
11321 vdw_lambda_2*tail_corr_dUdl_2);
11324 return scl*tail_corr_ener;
11334 return scl*(tail_corr_virial + vdw_lambda_1*tail_corr_virial_1 +
11335 vdw_lambda_2*tail_corr_virial_2);
11338 return scl*tail_corr_virial;
11343 #ifdef MEM_OPT_VERSION 11346 int Molecule::checkExclByIdx(
int idx1,
int atom1,
int atom2)
const {
11348 int amin = exclChkSigPool[idx1].min;
11349 int amax = exclChkSigPool[idx1].max;
11350 int dist21 = atom2 - atom1;
11351 if ( dist21 < amin || dist21 > amax )
return 0;
11352 else return exclChkSigPool[idx1].flags[dist21-amin];
11358 int amin = all_exclusions[atom1].min;
11359 int amax = all_exclusions[atom1].max;
11360 if ( atom2 < amin || atom2 > amax )
return 0;
11361 else return all_exclusions[atom1].flags[atom2-amin];
11381 is_lonepairs_psf = 0;
11389 read_parm(amber_data);
11391 #ifndef MEM_OPT_VERSION 11394 assignLCPOTypes( 1 );
11406 is_lonepairs_psf = 0;
11414 read_parm(amber_data);
11416 #ifndef MEM_OPT_VERSION 11419 assignLCPOTypes( 1 );
11424 #ifdef MEM_OPT_VERSION 11425 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
11427 int i,j,ntheth,nphih,current_index,a1,a2,
11428 max,min,index,found;
11431 NAMD_die(
"No data read from parm file yet!");
11434 numAtoms = amber_data->
Natom;
11435 atoms =
new Atom[numAtoms];
11447 if (atoms == NULL || atomNames == NULL )
11448 NAMD_die(
"memory allocation failed when reading atom information");
11451 for (i = 0; i < numAtoms; ++i) {
11452 const int resname_length = amber_data->
ResNames[amber_data->
AtomRes[i]].size();
11453 const int atomname_length = amber_data->
AtomNames[i].size();
11454 const int atomtype_length = amber_data->
AtomSym[i].size();
11455 atomNames[i].resname = nameArena->getNewArray(resname_length+1);
11456 atomNames[i].atomname = nameArena->getNewArray(atomname_length+1);
11457 atomNames[i].atomtype = nameArena->getNewArray(atomtype_length+1);
11458 if (atomNames[i].resname == NULL || atomNames[i].atomname == NULL || atomNames[i].atomtype == NULL)
11459 NAMD_die(
"memory allocation failed when reading atom information");
11461 strncpy(atomNames[i].resname, amber_data->
ResNames[amber_data->
AtomRes[i]].c_str(), resname_length+1);
11462 strncpy(atomNames[i].atomname, amber_data->
AtomNames[i].c_str(), atomname_length+1);
11463 strncpy(atomNames[i].atomtype, amber_data->
AtomSym[i].c_str(), atomtype_length+1);
11465 strtok(atomNames[i].resname,
" ");
11466 strtok(atomNames[i].atomname,
" ");
11467 strtok(atomNames[i].atomtype,
" ");
11468 atoms[i].mass = amber_data->
Masses[i];
11515 atoms[i].charge = amber_data->
Charges[i] / 18.2223;
11516 atoms[i].vdw_type = amber_data->
Iac[i] - 1;
11518 if ( tmpResLookup ) tmpResLookup =
11521 if(atomSegResids) {
11523 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
11530 }
else if (atoms[i].mass <= 0.05) {
11532 ++numZeroMassAtoms;
11533 }
else if (atoms[i].mass < 1.0) {
11539 }
else if (atoms[i].mass <=3.5) {
11541 }
else if ((atomNames[i].atomname[0] ==
'O') &&
11542 (atoms[i].mass >= 14.0) &&
11543 (atoms[i].mass <= 18.0)) {
11556 if (amber_data->
Nbonh + amber_data->
Nbona > 0) {
11558 if (bonds == NULL || amber_data->
Nbona < 0)
11559 NAMD_die(
"memory allocation failed when reading bond information");
11561 for (i=0; i<amber_data->
Nbonh; ++i)
11562 { bonds[numBonds].atom1 = amber_data->
BondHAt1[i] / 3;
11563 bonds[numBonds].atom2 = amber_data->
BondHAt2[i] / 3;
11564 bonds[numBonds].bond_type = amber_data->
BondHNum[i] - 1;
11565 if (bonds[numBonds].atom1>=numAtoms || bonds[numBonds].atom2>=numAtoms ||
11566 bonds[numBonds].bond_type>=amber_data->
Numbnd)
11567 {
char err_msg[128];
11568 sprintf(err_msg,
"BOND (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11576 for (i=amber_data->
Nbonh; i<amber_data->Nbonh+amber_data->
Nbona; ++i)
11577 { bonds[numBonds].atom1 = amber_data->
BondAt1[i-amber_data->
Nbonh] / 3;
11578 bonds[numBonds].atom2 = amber_data->
BondAt2[i-amber_data->
Nbonh] / 3;
11579 bonds[numBonds].bond_type = amber_data->
BondNum[i-amber_data->
Nbonh] - 1;
11580 if (bonds[i].atom1>=numAtoms || bonds[i].atom2>=numAtoms ||
11581 bonds[i].bond_type>=amber_data->
Numbnd)
11582 {
char err_msg[128];
11583 sprintf(err_msg,
"BOND (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-amber_data->
Nbonh);
11592 if ( numBonds != amber_data->
Nbonh + amber_data->
Nbona) {
11594 " bonds with zero force constants.\n" <<
endi;
11596 "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
endi;
11602 { ntheth = amber_data->
Ntheth;
11603 angles =
new Angle[numAngles];
11604 if (angles == NULL || numAngles < ntheth)
11605 NAMD_die(
"memory allocation failed when reading angle information");
11607 for (i=0; i<ntheth; ++i)
11608 { angles[i].atom1 = amber_data->
AngleHAt1[i] / 3;
11609 angles[i].atom2 = amber_data->
AngleHAt2[i] / 3;
11610 angles[i].atom3 = amber_data->
AngleHAt3[i] / 3;
11611 angles[i].angle_type = amber_data->
AngleHNum[i] - 1;
11612 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11613 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11614 {
char err_msg[128];
11615 sprintf(err_msg,
"ANGLE (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11620 for (i=ntheth; i<numAngles; ++i)
11621 { angles[i].atom1 = amber_data->
AngleAt1[i-ntheth] / 3;
11622 angles[i].atom2 = amber_data->
AngleAt2[i-ntheth] / 3;
11623 angles[i].atom3 = amber_data->
AngleAt3[i-ntheth] / 3;
11624 angles[i].angle_type = amber_data->
AngleNum[i-ntheth] - 1;
11625 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11626 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11627 {
char err_msg[128];
11628 sprintf(err_msg,
"ANGLE (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-ntheth);
11643 if (exclusions == NULL && amber_data->
Nnb > 0)
11644 NAMD_die(
"memory allocation failed when reading exclusion information");
11646 for (i=0; i<numAtoms; ++i)
11647 for (j=0; j<amber_data->
Iblo[i]; ++j)
11648 {
if (current_index >= amber_data->
Nnb)
11649 {
char err_msg[128];
11650 sprintf(err_msg,
"EXCLUSION INDEX EXCEEDS NUMBER OF EXLCUSIONS %d IN AMBER FILE, AT ATOM #%d\n",
11651 amber_data->
Nnb, i+1);
11656 if (amber_data->
ExclAt[current_index] != 0)
11659 a2 = amber_data->
ExclAt[current_index] - 1;
11665 sprintf(err_msg,
"Atom #%d has exclusion with atom #%d, in reverse order.", i+1, a2+1);
11669 {
char err_msg[128];
11670 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN AMBER FILE\n", i+1);
11673 else if (a2 >= numAtoms)
11674 {
char err_msg[128];
11675 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN AMBER FILE",
11676 a2+1, numAtoms, current_index+1);
11679 exclusions[numExclusions].atom1 = i;
11680 exclusions[numExclusions].atom2 = a2;
11685 if (current_index < amber_data->Nnb)
11686 {
char err_msg[128];
11687 sprintf(err_msg,
"Num of exclusions recorded (%d) is smaller than what it's supposed to be (%d)",
11688 current_index,amber_data->
Nnb);
11694 numDihedrals = amber_data->
Nphih + amber_data->
Nphia;
11695 if (numDihedrals > 0)
11696 { nphih = amber_data->
Nphih;
11697 dihedrals =
new Dihedral[numDihedrals];
11698 if (dihedrals == NULL || numDihedrals < nphih)
11699 NAMD_die(
"memory allocation failed when reading dihedral information");
11701 for (i=0; i<nphih; ++i)
11702 { dihedrals[i].atom1 = amber_data->
DihHAt1[i] / 3;
11703 dihedrals[i].atom2 = amber_data->
DihHAt2[i] / 3;
11704 dihedrals[i].atom3 = amber_data->
DihHAt3[i] / 3;
11705 dihedrals[i].atom4 = amber_data->
DihHAt4[i] / 3;
11706 dihedrals[i].dihedral_type = amber_data->
DihHNum[i] - 1;
11709 for (i=nphih; i<numDihedrals; ++i)
11710 { dihedrals[i].atom1 = amber_data->
DihAt1[i-nphih] / 3;
11711 dihedrals[i].atom2 = amber_data->
DihAt2[i-nphih] / 3;
11712 dihedrals[i].atom3 = amber_data->
DihAt3[i-nphih] / 3;
11713 dihedrals[i].atom4 = amber_data->
DihAt4[i-nphih] / 3;
11714 dihedrals[i].dihedral_type = amber_data->
DihNum[i-nphih] - 1;
11729 for (i=0; i<numDihedrals; ++i)
11730 {
if (dihedrals[i].atom3 < 0 || dihedrals[i].atom4 < 0)
11731 { dihedrals[i].atom3 = abs(dihedrals[i].atom3);
11732 dihedrals[i].atom4 = abs(dihedrals[i].atom4);
11735 {
if (dihedrals[i].atom1 < dihedrals[i].atom4)
11736 a1=dihedrals[i].atom1, a2=dihedrals[i].atom4;
11738 a1=dihedrals[i].atom4, a2=dihedrals[i].atom1;
11742 min=0, max=numExclusions-1;
11743 while (!found && min<=max)
11744 { index = (min+max)/2;
11745 if (exclusions[index].atom1 == a1)
11747 else if (exclusions[index].atom1 < a1)
11753 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11756 for (j=index-1; j>=0 && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; --j);
11757 if (j<0 || exclusions[j].atom1!=a1)
11758 for (j=index; j<numExclusions && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; ++j);
11759 if (j<numExclusions && exclusions[j].atom1==a1)
11760 exclusions[j].modified = 1;
11762 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11764 if (dihedrals[i].atom1>=numAtoms || dihedrals[i].atom2>=numAtoms ||
11765 dihedrals[i].atom3>=numAtoms || dihedrals[i].atom4>=numAtoms ||
11766 dihedrals[i].dihedral_type>=amber_data->
Nptra)
11767 {
char err_msg[128];
11768 sprintf(err_msg,
"DIHEDRAL # %d OVERFLOW IN PARM FILE", i+1);
11774 if (amber_data -> HasCMAP) {
11776 if (numCrossterms > 0) {
11777 crossterms =
new Crossterm[numCrossterms];
11779 for (i = 0; i < numCrossterms; ++i) {
11782 crossterms[i].atom2 = amber_data->
CMAPIndex[6*i+1]-1;
11783 crossterms[i].atom3 = amber_data->
CMAPIndex[6*i+2]-1;
11784 crossterms[i].atom4 = amber_data->
CMAPIndex[6*i+3]-1;
11786 crossterms[i].atom5 = amber_data->
CMAPIndex[6*i+1]-1;
11787 crossterms[i].atom6 = amber_data->
CMAPIndex[6*i+2]-1;
11788 crossterms[i].atom7 = amber_data->
CMAPIndex[6*i+3]-1;
11789 crossterms[i].atom8 = amber_data->
CMAPIndex[6*i+4]-1;
11791 crossterms[i].crossterm_type = amber_data->
CMAPIndex[6*i+5]-1;
11796 numRealBonds = numBonds;
11797 build_atom_status();
11813 void Molecule::read_parm(
Ambertoppar *amber_data)
11815 #ifdef MEM_OPT_VERSION 11816 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
11818 int i,j,ntheth,nphih,current_index,a1,a2,
11819 max,min,index,found;
11822 NAMD_die(
"No data read from parm file yet!");
11825 numAtoms = amber_data->
Natom;
11826 atoms =
new Atom[numAtoms];
11833 if (atoms == NULL || atomNames == NULL )
11834 NAMD_die(
"memory allocation failed when reading atom information");
11836 for (i=0; i<numAtoms; ++i)
11837 { atomNames[i].resname = nameArena->getNewArray(5);
11838 atomNames[i].atomname = nameArena->getNewArray(5);
11839 atomNames[i].atomtype = nameArena->getNewArray(5);
11840 if (atomNames[i].resname == NULL || atomNames[i].atomname == NULL || atomNames[i].atomtype == NULL)
11841 NAMD_die(
"memory allocation failed when reading atom information");
11842 for (j=0; j<4; ++j)
11843 { atomNames[i].resname[j] = amber_data->
ResNames[amber_data->
AtomRes[i]*4+j];
11844 atomNames[i].atomname[j] = amber_data->
AtomNames[i*4+j];
11845 atomNames[i].atomtype[j] = amber_data->
AtomSym[i*4+j];
11847 atomNames[i].resname[4] = atomNames[i].atomname[4] = atomNames[i].atomtype[4] =
'\0';
11848 strtok(atomNames[i].resname,
" ");
11849 strtok(atomNames[i].atomname,
" ");
11850 strtok(atomNames[i].atomtype,
" ");
11851 atoms[i].mass = amber_data->
Masses[i];
11853 atoms[i].charge = amber_data->
Charges[i] / 18.2223;
11854 atoms[i].vdw_type = amber_data->
Iac[i] - 1;
11857 if ( tmpResLookup ) tmpResLookup =
11860 if(atomSegResids) {
11862 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
11870 }
else if (atoms[i].mass <= 0.05) {
11871 ++numZeroMassAtoms;
11873 }
else if (atoms[i].mass < 1.0) {
11875 }
else if (atoms[i].mass <=3.5) {
11877 }
else if ((atomNames[i].atomname[0] ==
'O') &&
11878 (atoms[i].mass >= 14.0) &&
11879 (atoms[i].mass <= 18.0)) {
11892 if (amber_data->
Nbonh + amber_data->
Nbona > 0)
11894 if (bonds == NULL || amber_data->
Nbona < 0)
11895 NAMD_die(
"memory allocation failed when reading bond information");
11897 for (i=0; i<amber_data->
Nbonh; ++i)
11898 { bonds[numBonds].atom1 = amber_data->
BondHAt1[i] / 3;
11899 bonds[numBonds].atom2 = amber_data->
BondHAt2[i] / 3;
11900 bonds[numBonds].bond_type = amber_data->
BondHNum[i] - 1;
11901 if (bonds[numBonds].atom1>=numAtoms || bonds[numBonds].atom2>=numAtoms ||
11902 bonds[numBonds].bond_type>=amber_data->
Numbnd)
11903 {
char err_msg[128];
11904 sprintf(err_msg,
"BOND (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11912 for (i=amber_data->
Nbonh; i<amber_data->Nbonh+amber_data->
Nbona; ++i)
11913 { bonds[numBonds].atom1 = amber_data->
BondAt1[i-amber_data->
Nbonh] / 3;
11914 bonds[numBonds].atom2 = amber_data->
BondAt2[i-amber_data->
Nbonh] / 3;
11915 bonds[numBonds].bond_type = amber_data->
BondNum[i-amber_data->
Nbonh] - 1;
11916 if (bonds[i].atom1>=numAtoms || bonds[i].atom2>=numAtoms ||
11917 bonds[i].bond_type>=amber_data->
Numbnd)
11918 {
char err_msg[128];
11919 sprintf(err_msg,
"BOND (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-amber_data->
Nbonh);
11928 if ( numBonds != amber_data->
Nbonh + amber_data->
Nbona) {
11930 " bonds with zero force constants.\n" <<
endi;
11932 "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
endi;
11938 { ntheth = amber_data->
Ntheth;
11939 angles =
new Angle[numAngles];
11940 if (angles == NULL || numAngles < ntheth)
11941 NAMD_die(
"memory allocation failed when reading angle information");
11943 for (i=0; i<ntheth; ++i)
11944 { angles[i].atom1 = amber_data->
AngleHAt1[i] / 3;
11945 angles[i].atom2 = amber_data->
AngleHAt2[i] / 3;
11946 angles[i].atom3 = amber_data->
AngleHAt3[i] / 3;
11947 angles[i].angle_type = amber_data->
AngleHNum[i] - 1;
11948 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11949 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11950 {
char err_msg[128];
11951 sprintf(err_msg,
"ANGLE (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11956 for (i=ntheth; i<numAngles; ++i)
11957 { angles[i].atom1 = amber_data->
AngleAt1[i-ntheth] / 3;
11958 angles[i].atom2 = amber_data->
AngleAt2[i-ntheth] / 3;
11959 angles[i].atom3 = amber_data->
AngleAt3[i-ntheth] / 3;
11960 angles[i].angle_type = amber_data->
AngleNum[i-ntheth] - 1;
11961 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11962 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11963 {
char err_msg[128];
11964 sprintf(err_msg,
"ANGLE (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-ntheth);
11979 if (exclusions == NULL && amber_data->
Nnb > 0)
11980 NAMD_die(
"memory allocation failed when reading exclusion information");
11982 for (i=0; i<numAtoms; ++i)
11983 for (j=0; j<amber_data->
Iblo[i]; ++j)
11984 {
if (current_index >= amber_data->
Nnb)
11985 {
char err_msg[128];
11986 sprintf(err_msg,
"EXCLUSION INDEX EXCEEDS NUMBER OF EXLCUSIONS %d IN AMBER FILE, AT ATOM #%d\n",
11987 amber_data->
Nnb, i+1);
11992 if (amber_data->
ExclAt[current_index] != 0)
11995 a2 = amber_data->
ExclAt[current_index] - 1;
12001 sprintf(err_msg,
"Atom #%d has exclusion with atom #%d, in reverse order.", i+1, a2+1);
12005 {
char err_msg[128];
12006 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN AMBER FILE\n", i+1);
12009 else if (a2 >= numAtoms)
12010 {
char err_msg[128];
12011 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN AMBER FILE",
12012 a2+1, numAtoms, current_index+1);
12015 exclusions[numExclusions].atom1 = i;
12016 exclusions[numExclusions].atom2 = a2;
12021 if (current_index < amber_data->Nnb)
12022 {
char err_msg[128];
12023 sprintf(err_msg,
"Num of exclusions recorded (%d) is smaller than what it's supposed to be (%d)",
12024 current_index,amber_data->
Nnb);
12030 numDihedrals = amber_data->
Nphih + amber_data->
Nphia;
12031 if (numDihedrals > 0)
12032 { nphih = amber_data->
Nphih;
12033 dihedrals =
new Dihedral[numDihedrals];
12034 if (dihedrals == NULL || numDihedrals < nphih)
12035 NAMD_die(
"memory allocation failed when reading dihedral information");
12037 for (i=0; i<nphih; ++i)
12038 { dihedrals[i].atom1 = amber_data->
DihHAt1[i] / 3;
12039 dihedrals[i].atom2 = amber_data->
DihHAt2[i] / 3;
12040 dihedrals[i].atom3 = amber_data->
DihHAt3[i] / 3;
12041 dihedrals[i].atom4 = amber_data->
DihHAt4[i] / 3;
12042 dihedrals[i].dihedral_type = amber_data->
DihHNum[i] - 1;
12045 for (i=nphih; i<numDihedrals; ++i)
12046 { dihedrals[i].atom1 = amber_data->
DihAt1[i-nphih] / 3;
12047 dihedrals[i].atom2 = amber_data->
DihAt2[i-nphih] / 3;
12048 dihedrals[i].atom3 = amber_data->
DihAt3[i-nphih] / 3;
12049 dihedrals[i].atom4 = amber_data->
DihAt4[i-nphih] / 3;
12050 dihedrals[i].dihedral_type = amber_data->
DihNum[i-nphih] - 1;
12065 for (i=0; i<numDihedrals; ++i)
12066 {
if (dihedrals[i].atom3 < 0 || dihedrals[i].atom4 < 0)
12067 { dihedrals[i].atom3 = abs(dihedrals[i].atom3);
12068 dihedrals[i].atom4 = abs(dihedrals[i].atom4);
12071 {
if (dihedrals[i].atom1 < dihedrals[i].atom4)
12072 a1=dihedrals[i].atom1, a2=dihedrals[i].atom4;
12074 a1=dihedrals[i].atom4, a2=dihedrals[i].atom1;
12078 min=0, max=numExclusions-1;
12079 while (!found && min<=max)
12080 { index = (min+max)/2;
12081 if (exclusions[index].atom1 == a1)
12083 else if (exclusions[index].atom1 < a1)
12089 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
12092 for (j=index-1; j>=0 && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; --j);
12093 if (j<0 || exclusions[j].atom1!=a1)
12094 for (j=index; j<numExclusions && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; ++j);
12095 if (j<numExclusions && exclusions[j].atom1==a1)
12096 exclusions[j].modified = 1;
12098 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
12100 if (dihedrals[i].atom1>=numAtoms || dihedrals[i].atom2>=numAtoms ||
12101 dihedrals[i].atom3>=numAtoms || dihedrals[i].atom4>=numAtoms ||
12102 dihedrals[i].dihedral_type>=amber_data->
Nptra)
12103 {
char err_msg[128];
12104 sprintf(err_msg,
"DIHEDRAL # %d OVERFLOW IN PARM FILE", i+1);
12110 numRealBonds = numBonds;
12111 build_atom_status();
12130 read_parm(gromacsTopFile);
12132 #ifndef MEM_OPT_VERSION 12135 assignLCPOTypes( 3 );
12153 #ifdef MEM_OPT_VERSION 12154 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
12162 atoms =
new Atom[numAtoms];
12169 if (atoms == NULL || atomNames == NULL )
12170 NAMD_die(
"memory allocation failed when reading atom information");
12174 for (i=0; i<numAtoms; ++i) {
12175 char *resname = nameArena->getNewArray(11);
12176 char *atomname = nameArena->getNewArray(11);
12177 char *atomtype = nameArena->getNewArray(11);
12178 int resnum,typenum;
12181 if (resname == NULL || atomname == NULL || atomtype == NULL)
12182 NAMD_die(
"memory allocation failed when reading atom information");
12185 gf->
getAtom(i,&resnum,resname,atomname,atomtype,&typenum,
12188 atomNames[i].resname = resname;
12189 atomNames[i].atomname = atomname;
12190 atomNames[i].atomtype = atomtype;
12191 atoms[i].mass = mass;
12192 atoms[i].charge = charge;
12193 atoms[i].vdw_type = typenum;
12196 if ( tmpResLookup ) tmpResLookup =
12197 tmpResLookup->
append(
"MAIN", resnum+1, i);
12199 if(atomSegResids) {
12201 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
12202 one->
resid = resnum+1;
12211 }
else if (atoms[i].mass <= 0.05) {
12213 }
else if (atoms[i].mass < 1.0) {
12215 }
else if (atoms[i].mass <=3.5) {
12217 }
else if ((atomNames[i].atomname[0] ==
'O') &&
12218 (atoms[i].mass >= 14.0) &&
12219 (atoms[i].mass <= 18.0)) {
12226 bonds =
new Bond[numBonds];
12228 NAMD_die(
"memory allocation failed when reading bond information");
12229 for(i=0;i<numBonds;i++) {
12232 gf->
getBond(i,&atom1,&atom2,&type);
12233 bonds[i].atom1 = atom1;
12234 bonds[i].atom2 = atom2;
12235 bonds[i].bond_type = (
Index)type;
12240 angles =
new Angle[numAngles];
12241 if (angles == NULL)
12242 NAMD_die(
"memory allocation failed when reading angle information");
12243 for(i=0;i<numAngles;i++) {
12245 int atom1,atom2,atom3;
12246 gf->
getAngle(i,&atom1,&atom2,&atom3,&type);
12248 angles[i].atom1 = atom1;
12249 angles[i].atom2 = atom2;
12250 angles[i].atom3 = atom3;
12252 angles[i].angle_type=type;
12256 exclusions =
new Exclusion[numExclusions];
12320 dihedrals =
new Dihedral[numDihedrals];
12321 if (dihedrals == NULL)
12322 NAMD_die(
"memory allocation failed when reading dihedral information");
12323 for(i=0;i<numDihedrals;i++) {
12325 int atom1,atom2,atom3,atom4;
12326 gf->
getDihedral(i,&atom1,&atom2,&atom3,&atom4,&type);
12327 dihedrals[i].atom1 = atom1;
12328 dihedrals[i].atom2 = atom2;
12329 dihedrals[i].atom3 = atom3;
12330 dihedrals[i].atom4 = atom4;
12331 dihedrals[i].dihedral_type = type;
12344 const_cast<GromacsTopFile*
>(gf)->getPairLJArrays2(indxLJA, indxLJB, pairC6, pairC12);
12347 gromacsPair_type[i] = i;
12356 pointerToLJBeg =
new int[numAtoms];
12357 pointerToLJEnd =
new int[numAtoms];
12359 for(
int i=0; i < numAtoms; i++) {
12360 pointerToLJBeg[i] = -1;
12361 pointerToLJEnd[i] = -2;
12364 if(pointerToLJBeg[indxLJA[i]] == -1) {
12365 pointerToLJBeg[indxLJA[i]] = i;
12366 oldIndex = indxLJA[i];
12368 pointerToLJEnd[oldIndex] = i;
12381 const_cast<GromacsTopFile*
>(gf)->getPairGaussArrays2(indxGaussA, indxGaussB, gA, gMu1, giSigma1, gMu2, giSigma2, gRepulsive);
12384 pointerToGaussBeg =
new int[numAtoms];
12385 pointerToGaussEnd =
new int[numAtoms];
12386 for(
int i=0; i < numAtoms; i++) {
12387 pointerToGaussBeg[i] = -1;
12388 pointerToGaussEnd[i] = -2;
12392 if(pointerToGaussBeg[indxGaussA[i]] == -1) {
12393 pointerToGaussBeg[indxGaussA[i]] = i;
12394 oldIndex = indxGaussA[i];
12396 pointerToGaussEnd[oldIndex] = i;
12399 iout <<
iINFO <<
"Finished reading explicit pair from Gromacs file:\n" <<
12400 iINFO <<
"Found a total of: " << numPair <<
" explicit pairs--of which: " <<
12402 " are Gaussian style pairs.\n" <<
endi;
12406 #if GROMACS_EXCLUSIONS 12409 int* atom1 =
new int[numExclusions];
12410 int* atom2 =
new int[numExclusions];
12411 for(
int j=0; j<numExclusions;j++) {
12417 read_exclusions(atom1,atom2,numExclusions);
12481 numRealBonds = numBonds;
12482 build_atom_status();
12487 #ifndef MEM_OPT_VERSION 12501 #ifdef MEM_OPT_VERSION 12503 Index Molecule::insert_new_mass(
Real newMass){
12505 for(
int i=massPoolSize-1; i>=0; i--){
12506 if(fabs(atomMassPool[i]-newMass)<=1e-6)
12510 Real *tmp =
new Real[massPoolSize+1];
12511 tmp[massPoolSize] = newMass;
12512 memcpy((
void *)tmp, (
const void *)atomMassPool,
sizeof(
Real)*massPoolSize);
12513 delete [] atomMassPool;
12514 atomMassPool = tmp;
12516 return (
Index)(massPoolSize-1);
12519 void Molecule::addNewExclSigPool(
const vector<ExclusionSignature>& newExclSigPool){
12521 for(
int i=0; i<exclSigPoolSize; i++)
12522 tmpExclSigPool[i] = exclSigPool[i];
12523 for(
int i=0; i<newExclSigPool.size(); i++)
12524 tmpExclSigPool[i+exclSigPoolSize] = newExclSigPool[i];
12526 exclSigPoolSize += newExclSigPool.size();
12527 exclSigPool = tmpExclSigPool;
12531 msg->
put((
short)tupleType);
12532 msg->
put(numOffset);
12533 msg->
put(numOffset, offset);
12534 msg->
put(tupleParamType);
12543 msg->
get(numOffset);
12545 offset =
new int[numOffset];
12546 msg->
get(numOffset*
sizeof(
int), (
char *)offset);
12548 msg->
get(tupleParamType);
12554 for(
int i=0; i<bondCnt; i++)
12555 bondSigs[i].pack(msg);
12557 msg->
put(angleCnt);
12558 for(
int i=0; i<angleCnt; i++)
12559 angleSigs[i].pack(msg);
12561 msg->
put(dihedralCnt);
12562 for(
int i=0; i<dihedralCnt; i++)
12563 dihedralSigs[i].pack(msg);
12565 msg->
put(improperCnt);
12566 for(
int i=0; i<improperCnt; i++)
12567 improperSigs[i].pack(msg);
12569 msg->
put(crosstermCnt);
12570 for(
int i=0; i<crosstermCnt; i++)
12571 crosstermSigs[i].pack(msg);
12574 msg->
put(gromacsPairCnt);
12575 for(
int i=0; i<gromacsPairCnt; i++)
12576 gromacsPairSigs[i].pack(msg);
12581 delete [] bondSigs;
12584 for(
int i=0; i<bondCnt; i++)
12585 bondSigs[i].unpack(msg);
12586 }
else bondSigs = NULL;
12588 msg->
get(angleCnt);
12589 delete [] angleSigs;
12592 for(
int i=0; i<angleCnt; i++)
12593 angleSigs[i].unpack(msg);
12594 }
else angleSigs = NULL;
12596 msg->
get(dihedralCnt);
12597 delete [] dihedralSigs;
12600 for(
int i=0; i<dihedralCnt; i++)
12601 dihedralSigs[i].unpack(msg);
12602 }
else dihedralSigs = NULL;
12604 msg->
get(improperCnt);
12605 delete [] improperSigs;
12608 for(
int i=0; i<improperCnt; i++)
12609 improperSigs[i].unpack(msg);
12610 }
else improperSigs = NULL;
12612 msg->
get(crosstermCnt);
12613 delete [] crosstermSigs;
12614 if(crosstermCnt>0){
12616 for(
int i=0; i<crosstermCnt; i++)
12617 crosstermSigs[i].unpack(msg);
12618 }
else crosstermSigs = NULL;
12622 msg->
get(gromacsPairCnt);
12623 delete [] gromacsPairSigs;
12624 if(gromacsPairCnt>0){
12626 for(
int i=0; i<gromacsPairCnt; i++)
12627 gromacsPairSigs[i].unpack(msg);
12628 }
else gromacsPairSigs = NULL;
12642 origTupleCnt = bondCnt;
12643 tupleSigs= bondSigs;
12644 for(
int i=0; i<origTupleCnt; i++){
12649 delete [] tupleSigs;
12651 }
else if(bondCnt!=origTupleCnt){
12654 for(
int i=0; i<origTupleCnt; i++){
12656 newTupleSigs[idx] = tupleSigs[i];
12660 delete [] tupleSigs;
12661 bondSigs = newTupleSigs;
12667 origTupleCnt = angleCnt;
12668 tupleSigs = angleSigs;
12669 for(
int i=0; i<origTupleCnt; i++){
12674 delete [] tupleSigs;
12676 }
else if(angleCnt!=origTupleCnt){
12679 for(
int i=0; i<origTupleCnt; i++){
12681 newTupleSigs[idx] = tupleSigs[i];
12685 delete [] tupleSigs;
12686 angleSigs = newTupleSigs;
12692 origTupleCnt = dihedralCnt;
12693 tupleSigs = dihedralSigs;
12694 for(
int i=0; i<origTupleCnt; i++){
12698 if(dihedralCnt==0){
12699 delete [] tupleSigs;
12700 dihedralSigs = NULL;
12701 }
else if(dihedralCnt!=origTupleCnt){
12704 for(
int i=0; i<origTupleCnt; i++){
12706 newTupleSigs[idx] = tupleSigs[i];
12710 delete [] tupleSigs;
12711 dihedralSigs = newTupleSigs;
12718 origTupleCnt = improperCnt;
12719 tupleSigs = improperSigs;
12720 for(
int i=0; i<origTupleCnt; i++){
12724 if(improperCnt==0){
12725 delete [] tupleSigs;
12726 improperSigs = NULL;
12727 }
else if(improperCnt!=origTupleCnt){
12730 for(
int i=0; i<origTupleCnt; i++){
12732 newTupleSigs[idx] = tupleSigs[i];
12736 delete [] tupleSigs;
12737 improperSigs = newTupleSigs;
12743 origTupleCnt = crosstermCnt;
12744 tupleSigs = crosstermSigs;
12745 for(
int i=0; i<origTupleCnt; i++){
12749 if(crosstermCnt==0){
12750 delete [] tupleSigs;
12751 crosstermSigs = NULL;
12752 }
else if(crosstermCnt!=origTupleCnt){
12755 for(
int i=0; i<origTupleCnt; i++){
12757 newTupleSigs[idx] = tupleSigs[i];
12761 delete [] tupleSigs;
12762 crosstermSigs = newTupleSigs;
12769 origTupleCnt = gromacsPairCnt;
12770 tupleSigs = gromacsPairSigs;
12771 for(
int i=0; i<origTupleCnt; i++){
12775 if(gromacsPairCnt==0){
12776 delete [] tupleSigs;
12777 gromacsPairSigs = NULL;
12778 }
else if(gromacsPairCnt!=origTupleCnt){
12781 for(
int i=0; i<origTupleCnt; i++){
12783 newTupleSigs[idx] = tupleSigs[i];
12787 delete [] tupleSigs;
12788 gromacsPairSigs = newTupleSigs;
12798 for(
int i=0; i<fullExclCnt; i++){
12799 if(fullOffset[i]==0)
continue;
12804 delete [] fullOffset;
12806 }
else if(newCnt!=fullExclCnt){
12807 int *tmpOffset =
new int[newCnt];
12809 for(
int i=0; i<fullExclCnt; i++){
12810 if(fullOffset[i]==0)
continue;
12811 tmpOffset[newCnt] = fullOffset[i];
12814 delete [] fullOffset;
12815 fullOffset = tmpOffset;
12816 fullExclCnt = newCnt;
12821 for(
int i=0; i<modExclCnt; i++){
12822 if(modOffset[i]==0)
continue;
12827 delete [] modOffset;
12829 }
else if(newCnt!=modExclCnt){
12830 int *tmpOffset =
new int[newCnt];
12832 for(
int i=0; i<modExclCnt; i++){
12833 if(modOffset[i]==0)
continue;
12834 tmpOffset[newCnt] = modOffset[i];
12837 delete [] modOffset;
12838 modOffset = tmpOffset;
12839 modExclCnt = newCnt;
12854 int high = fullExclCnt-1;
12855 int mid = (low+high)/2;
12857 if(offset<fullOffset[mid]){
12859 mid = (high+low)/2;
12860 }
else if(offset>fullOffset[mid]){
12862 mid = (high+low)/2;
12868 if(retidx!=-1)
return retidx;
12872 high = modExclCnt-1;
12873 mid = (low+high)/2;
12875 if(offset<modOffset[mid]){
12877 mid = (high+low)/2;
12878 }
else if(offset>modOffset[mid]){
12880 mid = (high+low)/2;
12890 msg->
put(fullExclCnt);
12891 msg->
put(fullExclCnt, fullOffset);
12892 msg->
put(modExclCnt);
12893 msg->
put(modExclCnt, modOffset);
12897 msg->
get(fullExclCnt);
12898 delete [] fullOffset;
12899 fullOffset =
new int[fullExclCnt];
12900 msg->
get(fullExclCnt*
sizeof(
int), (
char *)fullOffset);
12901 msg->
get(modExclCnt);
12902 delete [] modOffset;
12903 modOffset =
new int[modExclCnt];
12904 msg->
get(modExclCnt*
sizeof(
int), (
char *)modOffset);
12905 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 12911 #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 *)