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;
210 lcpoParamType = NULL;
219 #ifdef MEM_OPT_VERSION 227 atomChargePool = NULL;
228 eachAtomCharge = NULL;
241 #ifndef MEM_OPT_VERSION 247 dihedralsByAtom=NULL;
248 impropersByAtom=NULL;
249 crosstermsByAtom=NULL;
251 gromacsPairByAtom=NULL;
259 #ifdef MEM_OPT_VERSION 261 exclChkSigPool = NULL;
263 eachAtomExclSig = NULL;
265 fixedAtomsSet = NULL;
266 constrainedAtomsSet = NULL;
268 exclusionsByAtom=NULL;
269 fullExclusionsByAtom=NULL;
270 modExclusionsByAtom=NULL;
277 #ifdef MEM_OPT_VERSION 284 exPressureAtomFlags=NULL;
285 rigidBondLengths=NULL;
300 consTorqueIndexes=NULL;
301 consTorqueParams=NULL;
302 consForceIndexes=NULL;
305 moleculeStartIndex=NULL;
310 alch_unpert_bonds = NULL;
311 alch_unpert_angles = NULL;
312 alch_unpert_dihedrals = NULL;
322 #ifndef MEM_OPT_VERSION 360 numExPressureAtoms=0;
362 numFixedRigidBonds=0;
363 numMultipleDihedrals=0;
364 numMultipleImpropers=0;
373 numCalcFullExclusions=0;
381 num_alch_unpert_Bonds = 0;
382 num_alch_unpert_Angles = 0;
383 num_alch_unpert_Dihedrals = 0;
457 #ifdef MEM_OPT_VERSION 459 read_mol_signatures(filename, param, cfgList);
461 read_psf_file(filename, param);
464 assignLCPOTypes( 0 );
478 #ifdef MEM_OPT_VERSION 479 NAMD_die(
"Sorry, plugin IO is not supported in the memory optimized version.");
483 int optflags = MOLFILE_BADOPTIONS;
484 molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*
sizeof(molfile_atom_t));
485 memset(atomarray, 0, natoms*
sizeof(molfile_atom_t));
488 int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
489 if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
491 NAMD_die(
"ERROR: plugin failed reading structure data");
493 if(optflags == MOLFILE_BADOPTIONS) {
495 NAMD_die(
"ERROR: plugin didn't initialize optional data flags");
497 if(optflags & MOLFILE_OCCUPANCY) {
498 setOccupancyData(atomarray);
500 if(optflags & MOLFILE_BFACTOR) {
501 setBFactorData(atomarray);
504 plgLoadAtomBasics(atomarray);
511 int *bondtype, nbondtypes;
513 if(pIOHdl->read_bonds!=NULL) {
514 if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder,
515 &bondtype, &nbondtypes, &bondtypename)){
516 NAMD_die(
"ERROR: failed reading bond information.");
521 plgLoadBonds(from,to);
525 int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
526 int ctermcols, ctermrows;
527 int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes;
528 int *impropertypes, numimpropertypes;
529 char **angletypenames, **dihedraltypenames, **impropertypenames;
531 plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
532 if(pIOHdl->read_angles!=NULL) {
533 if(pIOHdl->read_angles(pIOFileHdl,
534 &numAngles, &plgAngles,
535 &angletypes, &numangletypes, &angletypenames,
536 &numDihedrals, &plgDihedrals,
537 &dihedraltypes, &numdihedraltypes, &dihedraltypenames,
538 &numImpropers, &plgImpropers,
539 &impropertypes, &numimpropertypes, &impropertypenames,
540 &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) {
541 NAMD_die(
"ERROR: failed reading angle information.");
545 if(numAngles!=0) plgLoadAngles(plgAngles);
546 if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
547 if(numImpropers!=0) plgLoadImpropers(plgImpropers);
548 if(numCrossterms!=0) plgLoadCrossterms(plgCterms);
550 numRealBonds = numBonds;
554 assignLCPOTypes( 2 );
577 if (atomNames != NULL)
584 if (resLookup != NULL)
588 if (drudeConsts != NULL)
delete [] drudeConsts;
589 if (lphosts != NULL)
delete [] lphosts;
590 if (anisos != NULL)
delete [] anisos;
591 if (tholes != NULL)
delete [] tholes;
592 if (lphostIndexes != NULL)
delete [] lphostIndexes;
596 if (lcpoParamType != NULL)
delete [] lcpoParamType;
598 #ifdef MEM_OPT_VERSION 599 if(eachAtomSig)
delete [] eachAtomSig;
608 if (dihedrals != NULL)
611 if (impropers != NULL)
614 if (crossterms != NULL)
615 delete [] crossterms;
617 if (exclusions != NULL)
618 delete [] exclusions;
624 if (acceptors != NULL)
627 #ifdef MEM_OPT_VERSION 628 if(exclSigPool)
delete [] exclSigPool;
629 if(exclChkSigPool)
delete [] exclChkSigPool;
630 if(eachAtomExclSig)
delete [] eachAtomExclSig;
631 if(fixedAtomsSet)
delete fixedAtomsSet;
632 if(constrainedAtomsSet)
delete constrainedAtomsSet;
634 if (bondsByAtom != NULL)
635 delete [] bondsByAtom;
637 if (anglesByAtom != NULL)
638 delete [] anglesByAtom;
640 if (dihedralsByAtom != NULL)
641 delete [] dihedralsByAtom;
643 if (impropersByAtom != NULL)
644 delete [] impropersByAtom;
646 if (crosstermsByAtom != NULL)
647 delete [] crosstermsByAtom;
649 if (exclusionsByAtom != NULL)
650 delete [] exclusionsByAtom;
652 if (fullExclusionsByAtom != NULL)
653 delete [] fullExclusionsByAtom;
655 if (modExclusionsByAtom != NULL)
656 delete [] modExclusionsByAtom;
658 if (all_exclusions != NULL)
659 delete [] all_exclusions;
662 if (gromacsPairByAtom != NULL)
663 delete [] gromacsPairByAtom;
667 if (tholesByAtom != NULL)
668 delete [] tholesByAtom;
669 if (anisosByAtom != NULL)
670 delete [] anisosByAtom;
675 if (lcpoParamType != NULL)
676 delete [] lcpoParamType;
678 if (fixedAtomFlags != NULL)
679 delete [] fixedAtomFlags;
681 if (stirIndexes != NULL)
682 delete [] stirIndexes;
685 #ifdef MEM_OPT_VERSION 686 if(clusterSigs != NULL){
687 delete [] clusterSigs;
693 if (clusterSize != NULL)
694 delete [] clusterSize;
696 if (exPressureAtomFlags != NULL)
697 delete [] exPressureAtomFlags;
699 if (rigidBondLengths != NULL)
700 delete [] rigidBondLengths;
703 if (fepAtomFlags != NULL)
704 delete [] fepAtomFlags;
705 if (alch_unpert_bonds != NULL)
706 delete [] alch_unpert_bonds;
707 if (alch_unpert_angles != NULL)
708 delete [] alch_unpert_angles;
709 if (alch_unpert_dihedrals != NULL)
710 delete [] alch_unpert_dihedrals;
714 if (ssAtomFlags != NULL)
715 delete [] ssAtomFlags;
716 if (ss_vdw_type != NULL)
717 delete [] ss_vdw_type;
718 if (ss_index != NULL)
722 if (qmAtomGroup != NULL)
723 delete [] qmAtomGroup;
725 if (qmAtmIndx != NULL)
728 if (qmAtmChrg != NULL)
732 if (qmGrpNumBonds != NULL)
733 delete [] qmGrpNumBonds;
735 if (qmGrpSizes != NULL)
736 delete [] qmGrpSizes;
738 if (qmDummyBondVal != NULL)
739 delete [] qmDummyBondVal;
741 if (qmMMNumTargs != NULL)
742 delete [] qmMMNumTargs;
747 if (qmGrpChrg != NULL)
750 if (qmGrpMult != NULL)
753 if (qmMeMMindx != NULL)
754 delete [] qmMeMMindx;
756 if (qmMeQMGrp != NULL)
759 if (qmLSSSize != NULL)
762 if (qmLSSIdxs != NULL)
765 if (qmLSSMass != NULL)
768 if (qmLSSRefSize != NULL)
769 delete [] qmLSSRefSize;
771 if (qmLSSRefIDs != NULL)
772 delete [] qmLSSRefIDs;
774 if (qmLSSRefMass != NULL)
775 delete [] qmLSSRefMass;
777 if (qmMMBond != NULL) {
778 for(
int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
779 if (qmMMBond[grpIndx] != NULL)
780 delete [] qmMMBond[grpIndx];
785 if (qmGrpBonds != NULL) {
786 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
787 if (qmGrpBonds[grpIndx] != NULL)
788 delete [] qmGrpBonds[grpIndx];
790 delete [] qmGrpBonds;
793 if (qmMMBondedIndx != NULL) {
794 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
795 if (qmMMBondedIndx[grpIndx] != NULL)
796 delete [] qmMMBondedIndx[grpIndx];
798 delete [] qmMMBondedIndx;
801 if (qmMMChargeTarget != NULL) {
802 for(
int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
803 if (qmMMChargeTarget[grpIndx] != NULL)
804 delete [] qmMMChargeTarget[grpIndx];
806 delete [] qmMMChargeTarget;
809 if (qmElementArray != NULL){
810 for(
int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
811 if (qmElementArray[atmInd] != NULL)
812 delete [] qmElementArray[atmInd];
814 delete [] qmElementArray;
817 if (qmDummyElement != NULL){
818 for(
int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
819 if (qmDummyElement[atmInd] != NULL)
820 delete [] qmDummyElement[atmInd];
822 delete [] qmDummyElement;
825 if (qmCustPCSizes != NULL){
826 delete [] qmCustPCSizes;
829 if (qmCustomPCIdxs != NULL){
830 delete [] qmCustomPCIdxs;
833 if (cSMDindex != NULL) {
834 for(
int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
835 if (cSMDindex[grpIndx] != NULL)
836 delete [] cSMDindex[grpIndx];
841 if (cSMDpairs != NULL) {
842 for(
int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) {
843 if (cSMDpairs[instIndx] != NULL)
844 delete [] cSMDpairs[instIndx];
849 if (cSMDindxLen != NULL)
850 delete [] cSMDindxLen;
853 if (cSMDVels != NULL)
855 if (cSMDcoffs != NULL)
858 #ifndef MEM_OPT_VERSION 865 #ifndef MEM_OPT_VERSION 886 void Molecule::read_psf_file(
char *fname,
Parameters *params)
896 if ( (psf_file =
Fopen(fname,
"r")) == NULL)
898 sprintf(err_msg,
"UNABLE TO OPEN .psf FILE %s", fname);
914 sprintf(err_msg,
"EMPTY .psf FILE %s", fname);
922 sprintf(err_msg,
"UNABLE TO FIND \"PSF\" STRING IN PSF FILE %s",
931 iout <<
iWARN <<
"Reading PSF supporting DRUDE without " 932 "enabling the Drude model in the simulation config file\n" <<
endi;
951 sprintf(err_msg,
"MISSING EVERYTHING BUT PSF FROM %s", fname);
959 sprintf(err_msg,
"CAN NOT FIND \"NTITLE\" STRING IN PSF FILE %s",
964 sscanf(buffer,
"%d", &NumTitle);
980 sprintf(err_msg,
"FOUND EOF INSTEAD OF NATOM IN PSF FILE %s",
993 sprintf(err_msg,
"DIDN'T FIND \"NATOM\" IN PSF FILE %s",
999 sscanf(buffer,
"%d", &numAtoms);
1001 read_atoms(psf_file, params);
1014 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NBONDS IN PSF");
1020 NAMD_die(
"DID NOT FIND NBOND AFTER ATOM LIST IN PSF");
1024 sscanf(buffer,
"%d", &numBonds);
1027 read_bonds(psf_file);
1040 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NTHETA IN PSF");
1046 NAMD_die(
"DID NOT FIND NTHETA AFTER BOND LIST IN PSF");
1050 sscanf(buffer,
"%d", &numAngles);
1053 read_angles(psf_file, params);
1066 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NPHI IN PSF");
1072 NAMD_die(
"DID NOT FIND NPHI AFTER ANGLE LIST IN PSF");
1076 sscanf(buffer,
"%d", &numDihedrals);
1079 read_dihedrals(psf_file, params);
1092 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NIMPHI IN PSF");
1098 NAMD_die(
"DID NOT FIND NIMPHI AFTER ATOM LIST IN PSF");
1102 sscanf(buffer,
"%d", &numImpropers);
1105 read_impropers(psf_file, params);
1118 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NDON IN PSF");
1124 NAMD_die(
"DID NOT FIND NDON AFTER ATOM LIST IN PSF");
1128 sscanf(buffer,
"%d", &numDonors);
1131 read_donors(psf_file);
1144 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NACC IN PSF");
1150 NAMD_die(
"DID NOT FIND NACC AFTER ATOM LIST IN PSF");
1154 sscanf(buffer,
"%d", &numAcceptors);
1157 read_acceptors(psf_file);
1166 NAMD_die(
"EOF ENCOUNTERED LOOKING FOR NNB IN PSF FILE");
1171 sscanf(buffer,
"%d", &numExclusions);
1174 read_exclusions(psf_file);
1193 int is_found_numlp = 0;
1194 int is_found_numaniso = 0;
1195 int is_found_ncrterm = 0;
1196 while ( ! is_found ) {
1202 if (ret_code != 0) {
1205 else if ( (is_found_numlp =
NAMD_find_word(buffer,
"NUMLP")) > 0) {
1206 is_found = is_found_numlp;
1208 else if ( (is_found_numaniso =
NAMD_find_word(buffer,
"NUMANISO")) > 0) {
1209 is_found = is_found_numaniso;
1211 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1212 is_found = is_found_ncrterm;
1216 if (is_found_numlp) {
1219 sscanf(buffer,
"%d", &numLphosts);
1222 if (numLphosts == 0) {
1227 is_lonepairs_psf = 0;
1228 }
else if (
simParams->CUDASOAintegrate) {
1236 NAMD_die(
"FOUND LONE PAIR HOSTS IN PSF WITH \"LONEPAIRS\" DISABLED IN CONFIG FILE");
1240 if (numBonds) process_bonds(params);
1242 if (is_found_numlp) {
1244 if (numLphosts > 0) read_lphosts(psf_file);
1248 is_found_numaniso = 0;
1249 is_found_ncrterm = 0;
1250 while ( ! is_found ) {
1256 if (ret_code != 0) {
1259 else if ( (is_found_numaniso =
NAMD_find_word(buffer,
"NUMANISO")) > 0) {
1260 is_found = is_found_numaniso;
1262 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1263 is_found = is_found_ncrterm;
1268 if (is_found_numaniso && is_drude_psf) {
1271 sscanf(buffer,
"%d", &numAnisos);
1272 if (numAnisos > 0) read_anisos(psf_file);
1276 is_found_ncrterm = 0;
1277 while ( ! is_found ) {
1283 if (ret_code != 0) {
1286 else if ( (is_found_ncrterm =
NAMD_find_word(buffer,
"NCRTERM")) > 0) {
1287 is_found = is_found_ncrterm;
1291 else if (is_drude_psf ) {
1292 NAMD_die(
"DID NOT FIND REQUIRED NUMANISO IN DRUDE PSF FILE");
1294 else if (is_found_numaniso ) {
1295 NAMD_die(
"FOUND NUMANISO IN PSF FILE MISSING DRUDE DESIGNATION");
1298 if (is_found_ncrterm) {
1301 sscanf(buffer,
"%d", &numCrossterms);
1302 if (numCrossterms > 0) read_crossterms(psf_file, params);
1310 if (is_drude_psf && numDrudeAtoms) {
1312 Bond *newbonds =
new Bond[numBonds+numDrudeAtoms];
1313 memcpy(newbonds, bonds, numBonds*
sizeof(
Bond));
1316 int origNumBonds = numBonds;
1317 for (i=0; i < numAtoms; i++) {
1318 if (!is_drude(i))
continue;
1319 Bond *b = &(bonds[numBonds++]);
1323 atomNames[i-1].atomtype, atomNames[i].atomtype, b
1326 if (numBonds-origNumBonds != numDrudeAtoms) {
1327 NAMD_die(
"must have same number of Drude particles and parents");
1332 numRealBonds = numBonds;
1333 build_atom_status();
1354 void Molecule::read_atoms(FILE *fd,
Parameters *params)
1359 int last_atom_number=0;
1361 char segment_name[11];
1362 char residue_number[11];
1363 char residue_name[11];
1371 atoms =
new Atom[numAtoms];
1383 hydrogenGroup.resize(0);
1385 if (atoms == NULL || atomNames == NULL )
1387 NAMD_die(
"memory allocation failed in Molecule::read_atoms");
1393 while (atom_number < numAtoms)
1406 read_count=sscanf(buffer,
"%d %s %s %s %s %s %f %f",
1407 &atom_number, segment_name, residue_number,
1408 residue_name, atom_name, atom_type, &charge, &mass);
1409 if (mass <= 0.05) ++numZeroMassAtoms;
1412 if (read_count != 8)
1416 sprintf(err_msg,
"BAD ATOM LINE FORMAT IN PSF FILE IN ATOM LINE %d\nLINE=%s",
1417 last_atom_number+1, buffer);
1430 read_count=sscanf(buffer,
1433 "%*d %*s %*s %*s %*s %*s %*f %*f %*d %f %f", &alpha, &
thole);
1434 if (read_count != 2)
1438 sprintf(err_msg,
"BAD ATOM LINE FORMAT IN PSF FILE " 1439 "IN ATOM LINE %d\nLINE=%s", last_atom_number+1, buffer);
1442 drudeConsts[atom_number-1].alpha = alpha;
1443 drudeConsts[atom_number-1].thole =
thole;
1444 if (fabs(alpha) >= 1e-6) ++numDrudeAtoms;
1450 if ( sscanf(atom_type,
"%d", &atom_type_num) > 0 )
1452 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.");
1456 if (atom_number != last_atom_number+1)
1460 sprintf(err_msg,
"ATOM NUMBERS OUT OF ORDER AT ATOM #%d OF PSF FILE",
1461 last_atom_number+1);
1470 int reslength = strlen(residue_name)+1;
1471 int namelength = strlen(atom_name)+1;
1472 int typelength = strlen(atom_type)+1;
1474 atomNames[atom_number-1].resname = nameArena->getNewArray(reslength);
1475 atomNames[atom_number-1].atomname = nameArena->getNewArray(namelength);
1476 atomNames[atom_number-1].atomtype = nameArena->getNewArray(typelength);
1478 if (atomNames[atom_number-1].resname == NULL)
1480 NAMD_die(
"memory allocation failed in Molecule::read_atoms");
1484 strcpy(atomNames[atom_number-1].resname, residue_name);
1485 strcpy(atomNames[atom_number-1].atomname, atom_name);
1486 strcpy(atomNames[atom_number-1].atomtype, atom_type);
1487 atoms[atom_number-1].mass = mass;
1488 atoms[atom_number-1].charge = charge;
1492 if ( tmpResLookup ) tmpResLookup =
1493 tmpResLookup->
append(segment_name, atoi(residue_number), atom_number-1);
1497 memcpy(one->
segname, segment_name, strlen(segment_name)+1);
1498 one->
resid = atoi(residue_number);
1503 }
else if (atoms[atom_number-1].mass <= 0.05) {
1505 }
else if (atoms[atom_number-1].mass < 1.0) {
1506 atoms[atom_number-1].status |=
DrudeAtom;
1507 }
else if (atoms[atom_number-1].mass <=3.5) {
1509 }
else if ((atomNames[atom_number-1].atomname[0] ==
'O') &&
1510 (atoms[atom_number-1].mass >= 14.0) &&
1511 (atoms[atom_number-1].mass <= 18.0)) {
1517 &(atoms[atom_number-1]));
1536 void Molecule::read_bonds(FILE *fd)
1544 bonds=
new Bond[numBonds];
1548 NAMD_die(
"memory allocations failed in Molecule::read_bonds");
1552 while (num_read < numBonds)
1564 if (atom_nums[j] >= numAtoms)
1568 sprintf(err_msg,
"BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1574 Bond *b = &(bonds[num_read]);
1575 b->
atom1=atom_nums[0];
1576 b->
atom2=atom_nums[1];
1585 void Molecule::process_bonds(
Parameters *params) {
1588 int origNumBonds = numBonds;
1590 int numZeroFrcBonds = 0;
1592 int numDrudeBonds = 0;
1598 bool bond_found =
false;
1599 for (
int old_read=0; old_read < origNumBonds; ++old_read)
1602 Bond *b = &(bonds[num_read]);
1603 *b = bonds[old_read];
1604 atom_nums[0] = b->
atom1;
1605 atom_nums[1] = b->
atom2;
1610 atomNames[atom_nums[0]].atomtype,
1611 atomNames[atom_nums[1]].atomtype,
1625 if (!is_lp_bond && !bond_found) {
1627 snprintf(err_msg,
sizeof(err_msg),
1628 "UNABLE TO FIND BOND PARAMETERS FOR %s %s (ATOMS %i %i)",
1629 atomNames[atom_nums[0]].atomtype, atomNames[atom_nums[1]].atomtype,
1634 numZeroFrcBonds += (k == 0.);
1635 numLPBonds += is_lp_bond;
1636 numDrudeBonds += is_drude_bond;
1646 if ( is_lonepairs_psf ) {
1647 if (k == 0. || is_lp_bond) --numBonds;
1650 numLPBonds -= is_lp_bond;
1651 if (k == 0. && !is_lp_bond) --numBonds;
1660 if (k == 0. || is_lp_bond || is_drude_bond) --numBonds;
1666 if ( num_read != numBonds ) {
1667 NAMD_bug(
"num_read != numBonds in Molecule::process_bonds()");
1671 if ( numBonds != origNumBonds ) {
1672 if (numZeroFrcBonds) {
1673 iout <<
iWARN <<
"Ignored " << numZeroFrcBonds <<
1674 " bonds with zero force constants.\n" <<
1675 iWARN <<
"Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
1679 iout <<
iWARN <<
"Ignored " << numLPBonds <<
1680 " bonds with lone pairs.\n" <<
1681 iWARN <<
"Will infer lonepair bonds from LPhost entries.\n" <<
endi;
1683 if (numDrudeBonds) {
1684 iout <<
iWARN <<
"Ignored " << numDrudeBonds <<
1685 " bonds with Drude particles.\n" <<
1686 iWARN <<
"Will use polarizability to assign Drude bonds.\n" <<
endi;
1699 alch_unpert_bonds=
new Bond[num_alch_unpert_Bonds];
1701 if (alch_unpert_bonds == NULL) {
1702 NAMD_die(
"memory allocations failed in Molecule::read_alch_unpert_bonds");
1705 while (num_read < num_alch_unpert_Bonds) {
1706 for (j=0; j<2; j++) {
1709 if (atom_nums[j] >= numAtoms) {
1712 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);
1717 Bond *b = &(alch_unpert_bonds[num_read]);
1718 b->
atom1=atom_nums[0];
1719 b->
atom2=atom_nums[1];
1743 void Molecule::read_angles(FILE *fd,
Parameters *params)
1749 int origNumAngles = numAngles;
1751 angles=
new Angle[numAngles];
1755 NAMD_die(
"memory allocation failed in Molecule::read_angles");
1759 while (num_read < numAngles)
1772 if (atom_nums[j] >= numAtoms)
1776 sprintf(err_msg,
"ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1782 angles[num_read].atom1=atom_nums[0];
1783 angles[num_read].atom2=atom_nums[1];
1784 angles[num_read].atom3=atom_nums[2];
1789 atomNames[atom_nums[0]].atomtype,
1790 atomNames[atom_nums[1]].atomtype,
1791 atomNames[atom_nums[2]].atomtype,
1792 &(angles[num_read]),
simParams->alchOn ? -1 : 0);
1793 if ( angles[num_read].angle_type == -1 ) {
1794 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n" 1799 Real k, t0, k_ub, r_ub;
1800 if ( angles[num_read].angle_type == -1 ) { k = -1.; k_ub = -1.; }
else 1802 if ( k == 0. && k_ub == 0. ) --numAngles;
1807 if ( numAngles != origNumAngles ) {
1808 iout <<
iWARN <<
"Ignored " << origNumAngles - numAngles <<
1809 " angles with zero force constants.\n" <<
endi;
1822 alch_unpert_angles=
new Angle[num_alch_unpert_Angles];
1824 if (alch_unpert_angles == NULL) {
1825 NAMD_die(
"memory allocation failed in Molecule::read_alch_unpert_angles");
1828 while (num_read < num_alch_unpert_Angles) {
1829 for (j=0; j<3; j++) {
1832 if (atom_nums[j] >= numAtoms) {
1834 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);
1839 alch_unpert_angles[num_read].atom1=atom_nums[0];
1840 alch_unpert_angles[num_read].atom2=atom_nums[1];
1841 alch_unpert_angles[num_read].atom3=atom_nums[2];
1864 void Molecule::read_dihedrals(FILE *fd,
Parameters *params)
1867 int last_atom_nums[4];
1871 Bool duplicate_bond;
1876 last_atom_nums[j] = -1;
1879 dihedrals =
new Dihedral[numDihedrals];
1881 if (dihedrals == NULL)
1883 NAMD_die(
"memory allocation failed in Molecule::read_dihedrals");
1887 while (num_read < numDihedrals)
1889 duplicate_bond =
TRUE;
1901 if (atom_nums[j] >= numAtoms)
1905 sprintf(err_msg,
"DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1910 if (atom_nums[j] != last_atom_nums[j])
1912 duplicate_bond =
FALSE;
1915 last_atom_nums[j] = atom_nums[j];
1925 if (multiplicity == 2)
1927 numMultipleDihedrals++;
1937 dihedrals[num_unique-1].atom1=atom_nums[0];
1938 dihedrals[num_unique-1].atom2=atom_nums[1];
1939 dihedrals[num_unique-1].atom3=atom_nums[2];
1940 dihedrals[num_unique-1].atom4=atom_nums[3];
1944 atomNames[atom_nums[0]].atomtype,
1945 atomNames[atom_nums[1]].atomtype,
1946 atomNames[atom_nums[2]].atomtype,
1947 atomNames[atom_nums[3]].atomtype,
1948 &(dihedrals[num_unique-1]),
1949 multiplicity,
simParams->alchOn ? -1 : 0);
1950 if ( dihedrals[num_unique-1].dihedral_type == -1 ) {
1951 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n" 1958 numDihedrals = num_unique;
1970 alch_unpert_dihedrals =
new Dihedral[num_alch_unpert_Dihedrals];
1972 if (alch_unpert_dihedrals == NULL) {
1973 NAMD_die(
"memory allocation failed in Molecule::read_alch_unpert_dihedrals");
1976 while (num_read < num_alch_unpert_Dihedrals) {
1977 for (j=0; j<4; j++) {
1980 if (atom_nums[j] >= numAtoms) {
1983 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);
1988 alch_unpert_dihedrals[num_read].atom1=atom_nums[0];
1989 alch_unpert_dihedrals[num_read].atom2=atom_nums[1];
1990 alch_unpert_dihedrals[num_read].atom3=atom_nums[2];
1991 alch_unpert_dihedrals[num_read].atom4=atom_nums[3];
2014 void Molecule::read_impropers(FILE *fd,
Parameters *params)
2017 int last_atom_nums[4];
2021 Bool duplicate_bond;
2027 last_atom_nums[j] = -1;
2030 impropers=
new Improper[numImpropers];
2032 if (impropers == NULL)
2034 NAMD_die(
"memory allocation failed in Molecule::read_impropers");
2038 while (num_read < numImpropers)
2040 duplicate_bond =
TRUE;
2052 if (atom_nums[j] >= numAtoms)
2056 sprintf(err_msg,
"IMPROPERS INDEX %d GREATER THAN NATOM %d IN IMPROPERS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
2060 if (atom_nums[j] != last_atom_nums[j])
2062 duplicate_bond =
FALSE;
2065 last_atom_nums[j] = atom_nums[j];
2076 if (multiplicity == 2)
2079 numMultipleImpropers++;
2090 impropers[num_unique-1].atom1=atom_nums[0];
2091 impropers[num_unique-1].atom2=atom_nums[1];
2092 impropers[num_unique-1].atom3=atom_nums[2];
2093 impropers[num_unique-1].atom4=atom_nums[3];
2097 atomNames[atom_nums[0]].atomtype,
2098 atomNames[atom_nums[1]].atomtype,
2099 atomNames[atom_nums[2]].atomtype,
2100 atomNames[atom_nums[3]].atomtype,
2101 &(impropers[num_unique-1]),
2110 numImpropers = num_unique;
2130 void Molecule::read_crossterms(FILE *fd,
Parameters *params)
2134 int last_atom_nums[8];
2137 Bool duplicate_bond;
2142 last_atom_nums[j] = -1;
2145 crossterms=
new Crossterm[numCrossterms];
2147 if (crossterms == NULL)
2149 NAMD_die(
"memory allocation failed in Molecule::read_crossterms");
2153 while (num_read < numCrossterms)
2155 duplicate_bond =
TRUE;
2167 if (atom_nums[j] >= numAtoms)
2171 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);
2175 if (atom_nums[j] != last_atom_nums[j])
2177 duplicate_bond =
FALSE;
2180 last_atom_nums[j] = atom_nums[j];
2186 iout <<
iWARN <<
"Duplicate cross-term detected.\n" <<
endi;
2190 crossterms[num_read].atom1=atom_nums[0];
2191 crossterms[num_read].atom2=atom_nums[1];
2192 crossterms[num_read].atom3=atom_nums[2];
2193 crossterms[num_read].atom4=atom_nums[3];
2194 crossterms[num_read].atom5=atom_nums[4];
2195 crossterms[num_read].atom6=atom_nums[5];
2196 crossterms[num_read].atom7=atom_nums[6];
2197 crossterms[num_read].atom8=atom_nums[7];
2201 atomNames[atom_nums[0]].atomtype,
2202 atomNames[atom_nums[1]].atomtype,
2203 atomNames[atom_nums[2]].atomtype,
2204 atomNames[atom_nums[3]].atomtype,
2205 atomNames[atom_nums[4]].atomtype,
2206 atomNames[atom_nums[5]].atomtype,
2207 atomNames[atom_nums[6]].atomtype,
2208 atomNames[atom_nums[7]].atomtype,
2209 &(crossterms[num_read]));
2211 if(!duplicate_bond) num_read++;
2214 numCrossterms = num_read;
2237 void Molecule::read_donors(FILE *fd)
2245 donors=
new Bond[numDonors];
2249 NAMD_die(
"memory allocations failed in Molecule::read_donors");
2253 while (num_read < numDonors)
2265 if (d[j] >= numAtoms)
2270 "DONOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE",
2271 d[j]+1, numAtoms, num_read+1);
2281 Bond *b = &(donors[num_read]);
2311 void Molecule::read_acceptors(FILE *fd)
2319 acceptors=
new Bond[numAcceptors];
2321 if (acceptors == NULL)
2323 NAMD_die(
"memory allocations failed in Molecule::read_acceptors");
2327 while (num_read < numAcceptors)
2339 if (d[j] >= numAtoms)
2343 sprintf(err_msg,
"ACCEPTOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE", d[j]+1, numAtoms, num_read+1);
2353 Bond *b = &(acceptors[num_read]);
2392 void Molecule::read_exclusions(FILE *fd)
2395 int *exclusion_atoms;
2396 register int num_read=0;
2399 register int insert_index=0;
2403 exclusions =
new Exclusion[numExclusions];
2404 exclusion_atoms =
new int[numExclusions];
2406 if ( (exclusions == NULL) || (exclusion_atoms == NULL) )
2408 NAMD_die(
"memory allocation failed in Molecule::read_exclusions");
2412 for (num_read=0; num_read<numExclusions; num_read++)
2420 if (exclusion_atoms[num_read] >= numAtoms)
2424 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);
2433 for (num_read=0; num_read<numAtoms; num_read++)
2439 if (current_index>numExclusions)
2443 sprintf(err_msg,
"EXCLUSION INDEX %d LARGER THAN NUMBER OF EXLCUSIONS %d IN PSF FILE, EXCLUSION #%d\n",
2444 current_index+1, numExclusions, num_read);
2451 if (current_index != last_index)
2457 for (insert_index=last_index;
2458 insert_index<current_index; insert_index++)
2465 int a2 = exclusion_atoms[insert_index];
2467 exclusions[insert_index].atom1 = a1;
2468 exclusions[insert_index].atom2 = a2;
2469 }
else if ( a2 < a1 ) {
2470 exclusions[insert_index].atom1 = a2;
2471 exclusions[insert_index].atom2 = a1;
2474 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN PSF FILE\n", a1+1);
2479 last_index=current_index;
2484 delete [] exclusion_atoms;
2500 void Molecule::read_exclusions(
int* atom_i,
int* atom_j,
int num_exclusion)
2504 exclusions =
new Exclusion[num_exclusion];
2505 int loop_counter = 0;
2509 if ( (exclusions == NULL) )
2511 NAMD_die(
"memory allocation failed in Molecule::read_exclusions");
2515 for (loop_counter = 0; loop_counter < num_exclusion; loop_counter++) {
2517 if ( (atom_i == NULL) || (atom_j == NULL) ) {
2518 NAMD_die(
"null pointer expection in Molecule::read_exclusions");
2521 a = atom_i[loop_counter];
2522 b = atom_j[loop_counter];
2524 exclusions[loop_counter].atom1 = a;
2525 exclusions[loop_counter].atom2 = b;
2527 exclusions[loop_counter].atom1 = b;
2528 exclusions[loop_counter].atom2 = a;
2530 exclusionSet.add(
Exclusion(exclusions[loop_counter].atom1,exclusions[loop_counter].atom2));
2534 iout <<
iINFO <<
"ADDED " << num_exclusion <<
" EXPLICIT EXCLUSIONS: THIS VALUE WILL *NOT* BE ADDED TO THE STRUCTURE SUMMARY\n" <<
endi;
2558 void Molecule::read_lphosts(FILE *fd)
2570 Real value1, value2, value3;
2573 lphosts =
new Lphost[numLphosts];
2574 if (lphosts == NULL) {
2575 NAMD_die(
"memory allocation failed in Molecule::read_lphosts");
2577 for (i = 0; i < numLphosts; i++) {
2580 read_count=sscanf(buffer,
"%d %d %6s %f %f %f",
2581 &numhosts, &index, weight, &value1, &value2, &value3);
2586 if (read_count != 6 || index != next_index ||
2587 strcmp(weight,
"F") != 0 || numhosts < 2 || 3 < numhosts) {
2589 sprintf(err_msg,
"BAD FORMAT FOR LPHOST LINE %d IN PSF FILE LINE\n" 2590 "LINE=%s\n", i+1, buffer);
2593 lphosts[i].numhosts = numhosts;
2594 next_index += numhosts + 1;
2595 if (numhosts == 2) {
2596 lphosts[i].distance = value1;
2597 lphosts[i].angle = value2;
2598 lphosts[i].dihedral = 0.0;
2601 lphosts[i].distance = value1;
2602 lphosts[i].angle = value2 * (
M_PI/180);
2603 lphosts[i].dihedral = value3 * (
M_PI/180);
2607 Bond *newbonds =
new Bond[numBonds+numLphosts];
2608 memcpy(newbonds, bonds, numBonds*
sizeof(
Bond));
2611 for (i = 0; i < numLphosts; i++) {
2617 lphosts[i].atom4 = ( lphosts[i].numhosts == 3 ?
2620 Bond *b = &(bonds[numBonds++]);
2621 b->
atom1 = lphosts[i].atom1;
2622 b->
atom2 = lphosts[i].atom2;
2637 void Molecule::read_anisos(FILE *fd)
2640 int numhosts, index, i, read_count;
2643 anisos =
new Aniso[numAnisos];
2646 NAMD_die(
"memory allocation failed in Molecule::read_anisos");
2648 for (i = 0; i < numAnisos; i++)
2652 read_count=sscanf(buffer,
"%f %f %f", &k11, &k22, &k33);
2653 if (read_count != 3)
2656 sprintf(err_msg,
"BAD FORMAT FOR ANISO LINE %d IN PSF FILE LINE\n" 2657 "LINE=%s\n", i+1, buffer);
2660 anisos[i].k11 = k11;
2661 anisos[i].k22 = k22;
2662 anisos[i].k33 = k33;
2664 for (i = 0; i < numAnisos; i++) {
2677 if (atomType[0] ==
'H' || atomType[0] ==
'h') {
2681 }
else if (atomType[0] ==
'C' || atomType[0] ==
'c') {
2684 strcmp(atomType,
"CT" )==0 )
2692 else if (numBonds == 2)
2694 else if (numBonds == 3)
2696 else if (numBonds == 4)
2704 else if (numBonds == 3)
2711 }
else if (atomType[0] ==
'N' || atomType[0] ==
'n') {
2712 if ( strcmp(atomType,
"N3" ) == 0 ) {
2715 else if (numBonds == 2)
2717 else if (numBonds == 3)
2725 else if (numBonds == 2)
2727 else if (numBonds == 3)
2734 }
else if (atomType[0] ==
'O' || atomType[0] ==
'o') {
2736 if ( strcmp(atomType,
"O" )==0) {
2738 }
else if (strcmp(atomType,
"O2" )==0) {
2743 else if (numBonds == 2)
2750 }
else if (atomType[0] ==
'S' || atomType[0] ==
's') {
2751 if ( strcmp(atomType,
"SH" )==0) {
2758 }
else if (atomType[0] ==
'P' || atomType[0] ==
'p') {
2761 else if (numBonds == 4)
2765 }
else if (atomType[0] ==
'Z') {
2767 }
else if ( strcmp(atomType,
"MG" )==0) {
2778 if (atomType[0] ==
'H') {
2782 }
else if (atomType[0] ==
'C') {
2784 atomType[1] ==
'T' ||
2785 strcmp(atomType,
"CP1" )==0 ||
2786 strcmp(atomType,
"CP2" )==0 ||
2787 strcmp(atomType,
"CP3" )==0 ||
2788 strcmp(atomType,
"CS" )==0 ) {
2791 else if (numBonds == 2)
2793 else if (numBonds == 3)
2795 else if (numBonds == 4)
2801 strcmp(atomType,
"C" )==0 ||
2802 strcmp(atomType,
"CA" )==0 ||
2803 strcmp(atomType,
"CC" )==0 ||
2804 strcmp(atomType,
"CD" )==0 ||
2805 strcmp(atomType,
"CN" )==0 ||
2806 strcmp(atomType,
"CY" )==0 ||
2807 strcmp(atomType,
"C3" )==0 ||
2808 strcmp(atomType,
"CE1" )==0 ||
2809 strcmp(atomType,
"CE2" )==0 ||
2810 strcmp(atomType,
"CST" )==0 ||
2811 strcmp(atomType,
"CAP" )==0 ||
2812 strcmp(atomType,
"COA" )==0 ||
2813 strcmp(atomType,
"CPT" )==0 ||
2814 strcmp(atomType,
"CPH1")==0 ||
2815 strcmp(atomType,
"CPH2")==0
2819 else if (numBonds == 3)
2828 }
else if (atomType[0] ==
'N') {
2833 strcmp(atomType,
"NH3" )==0 ||
2836 strcmp(atomType,
"NP" )==0
2840 else if (numBonds == 2)
2842 else if (numBonds == 3)
2848 strcmp(atomType,
"NY" )==0 ||
2849 strcmp(atomType,
"NC2" )==0 ||
2850 strcmp(atomType,
"N" )==0 ||
2851 strcmp(atomType,
"NH1" )==0 ||
2852 strcmp(atomType,
"NH2" )==0 ||
2853 strcmp(atomType,
"NR1" )==0 ||
2854 strcmp(atomType,
"NR2" )==0 ||
2855 strcmp(atomType,
"NR3" )==0 ||
2856 strcmp(atomType,
"NPH" )==0 ||
2857 strcmp(atomType,
"NC" )==0
2861 else if (numBonds == 2)
2863 else if (numBonds == 3)
2872 }
else if (atomType[0] ==
'O') {
2874 strcmp(atomType,
"OH1" )==0 ||
2875 strcmp(atomType,
"OS" )==0 ||
2876 strcmp(atomType,
"OC" )==0 ||
2877 strcmp(atomType,
"OT" )==0
2881 else if (numBonds == 2)
2886 strcmp(atomType,
"O" )==0 ||
2887 strcmp(atomType,
"OB" )==0 ||
2888 strcmp(atomType,
"OST" )==0 ||
2889 strcmp(atomType,
"OCA" )==0 ||
2890 strcmp(atomType,
"OM" )==0
2894 strcmp(atomType,
"OC" )==0
2902 }
else if (atomType[0] ==
'S') {
2909 }
else if (atomType[0] ==
'P') {
2912 else if (numBonds == 4)
2927 void Molecule::assignLCPOTypes(
int inputType) {
2928 int *heavyBonds =
new int[numAtoms];
2929 for (
int i = 0; i < numAtoms; i++)
2931 for (
int i = 0; i < numBonds; i++ ) {
2932 Bond curBond = bonds[i];
2933 int a1 = bonds[i].
atom1;
2934 int a2 = bonds[i].atom2;
2935 if (atoms[a1].mass > 2.f && atoms[a2].mass > 2.f) {
2941 lcpoParamType =
new int[numAtoms];
2944 for (
int i = 0; i < numAtoms; i++) {
2947 if (inputType == 1) {
2961 if ( atoms[i].mass < 1.5 && lcpoParamType[i] != 0 ) {
2964 lcpoParamType[i] = 0;
2968 lcpoParamType[i] = 0;
2971 CkPrintf(
"ERROR in Molecule::assignLCPOTypes(): " 2972 "Light atom given heavy atom LCPO type.\n");
2980 iout <<
iWARN <<
"LONE PAIRS TO BE IGNORED BY SASA\n" <<
endi;
2983 iout <<
iWARN <<
"DRUDE PARTICLES TO BE IGNORED BY SASA\n" <<
endi;
2986 delete [] heavyBonds;
2990 void Molecule::plgLoadAtomBasics(molfile_atom_t *atomarray){
2991 atoms =
new Atom[numAtoms];
2996 hydrogenGroup.resize(0);
3000 for(
int i=0; i<numAtoms; i++) {
3001 int reslength = strlen(atomarray[i].resname)+1;
3002 int namelength = strlen(atomarray[i].name)+1;
3003 int typelength = strlen(atomarray[i].type)+1;
3004 atomNames[i].resname = nameArena->getNewArray(reslength);
3005 atomNames[i].atomname = nameArena->getNewArray(namelength);
3006 atomNames[i].atomtype = nameArena->getNewArray(typelength);
3007 strcpy(atomNames[i].resname, atomarray[i].resname);
3008 strcpy(atomNames[i].atomname, atomarray[i].name);
3009 strcpy(atomNames[i].atomtype, atomarray[i].type);
3011 atoms[i].mass = atomarray[i].mass;
3012 atoms[i].charge = atomarray[i].charge;
3017 tmpResLookup = tmpResLookup->
append(atomarray[i].segid, atomarray[i].resid, i);
3022 memcpy(one->
segname, atomarray[i].segid, strlen(atomarray[i].segid)+1);
3023 one->
resid = atomarray[i].resid;
3027 }
else if(atoms[i].mass <= 0.05) {
3029 }
else if(atoms[i].mass < 1.0) {
3031 }
else if(atoms[i].mass <= 3.5) {
3033 }
else if((atomNames[i].atomname[0] ==
'O') &&
3034 (atoms[i].mass>=14.0) && (atoms[i].mass<=18.0)){
3042 void Molecule::plgLoadBonds(
int *from,
int *to){
3043 bonds =
new Bond[numBonds];
3044 int realNumBonds = 0;
3045 for(
int i=0; i<numBonds; i++) {
3046 Bond *thisBond = bonds+realNumBonds;
3047 thisBond->
atom1 = from[i]-1;
3048 thisBond->
atom2 = to[i]-1;
3051 atomNames[thisBond->
atom1].atomtype,
3052 atomNames[thisBond->
atom2].atomtype,
3058 if (is_lonepairs_psf) {
3060 if(k!=0. || is_lp(thisBond->
atom1) ||
3061 is_lp(thisBond->
atom2)) {
3065 if(k != 0.) realNumBonds++;
3069 if(numBonds != realNumBonds) {
3070 iout <<
iWARN <<
"Ignored" << numBonds-realNumBonds <<
3071 "bonds with zero force constants.\n" <<
endi;
3072 iout <<
iWARN <<
"Will get H-H distance in rigid H20 from H-O-H angle.\n" <<
endi;
3074 numBonds = realNumBonds;
3077 void Molecule::plgLoadAngles(
int *plgAngles)
3079 angles=
new Angle[numAngles];
3080 int *atomid = plgAngles;
3081 int numRealAngles = 0;
3082 for(
int i=0; i<numAngles; i++) {
3083 Angle *thisAngle = angles+numRealAngles;
3084 thisAngle->
atom1 = atomid[0]-1;
3085 thisAngle->
atom2 = atomid[1]-1;
3086 thisAngle->
atom3 = atomid[2]-1;
3090 atomNames[thisAngle->
atom1].atomtype,
3091 atomNames[thisAngle->
atom2].atomtype,
3092 atomNames[thisAngle->
atom3].atomtype,
3095 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n" 3099 Real k, t0, k_ub, r_ub;
3100 if ( thisAngle->
angle_type == -1 ) { k = -1.; k_ub = -1.; }
else 3102 if(k!=0. || k_ub!=0.) numRealAngles++;
3105 if(numAngles != numRealAngles) {
3106 iout <<
iWARN <<
"Ignored" << numAngles-numRealAngles <<
3107 " angles with zero force constants.\n" <<
endi;
3109 numAngles = numRealAngles;
3112 void Molecule::plgLoadDihedrals(
int *plgDihedrals)
3114 std::map< std::string, int > cache;
3117 int multiplicity = 1;
3119 lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
3120 dihedrals =
new Dihedral[numDihedrals];
3121 int numRealDihedrals = 0;
3122 int *atomid = plgDihedrals;
3123 for(
int i=0; i<numDihedrals; i++, atomid+=4) {
3124 Dihedral *thisDihedral = dihedrals + numRealDihedrals;
3126 for(
int j=0; j<4; j++) {
3127 if(atomid[j] != lastAtomIds[j]) {
3128 duplicate_bond =
FALSE;
3130 lastAtomIds[j] = atomid[j];
3133 if(duplicate_bond) {
3135 if(multiplicity==2) {
3136 numMultipleDihedrals++;
3143 thisDihedral->
atom1 = atomid[0]-1;
3144 thisDihedral->
atom2 = atomid[1]-1;
3145 thisDihedral->
atom3 = atomid[2]-1;
3146 thisDihedral->
atom4 = atomid[3]-1;
3149 sprintf(query,
"%10s :: %10s :: %10s :: %10s :: %d",
3150 atomNames[atomid[0]-1].atomtype,
3151 atomNames[atomid[1]-1].atomtype,
3152 atomNames[atomid[2]-1].atomtype,
3153 atomNames[atomid[3]-1].atomtype,
3155 auto search = cache.find(query);
3156 if ( search != cache.end() ) {
3160 atomNames[atomid[0]-1].atomtype,
3161 atomNames[atomid[1]-1].atomtype,
3162 atomNames[atomid[2]-1].atomtype,
3163 atomNames[atomid[3]-1].atomtype,
3164 thisDihedral, multiplicity,
simParams->alchOn ? -1 : 0);
3166 iout <<
iWARN <<
"ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n" 3173 numDihedrals = numRealDihedrals;
3176 void Molecule::plgLoadImpropers(
int *plgImpropers)
3179 int multiplicity = 1;
3181 lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
3182 impropers =
new Improper[numImpropers];
3183 int numRealImpropers = 0;
3184 int *atomid = plgImpropers;
3185 for(
int i=0; i<numImpropers; i++, atomid+=4) {
3186 Improper *thisImproper = impropers + numRealImpropers;
3188 for(
int j=0; j<4; j++) {
3189 if(atomid[j] != lastAtomIds[j]) {
3190 duplicate_bond =
FALSE;
3192 lastAtomIds[j] = atomid[j];
3195 if(duplicate_bond) {
3197 if(multiplicity==2) {
3198 numMultipleImpropers++;
3205 thisImproper->
atom1 = atomid[0]-1;
3206 thisImproper->
atom2 = atomid[1]-1;
3207 thisImproper->
atom3 = atomid[2]-1;
3208 thisImproper->
atom4 = atomid[3]-1;
3211 atomNames[atomid[0]-1].atomtype,
3212 atomNames[atomid[1]-1].atomtype,
3213 atomNames[atomid[2]-1].atomtype,
3214 atomNames[atomid[3]-1].atomtype,
3215 thisImproper, multiplicity);
3218 numImpropers = numRealImpropers;
3221 void Molecule::plgLoadCrossterms(
int *plgCterms)
3225 for(
int i=0; i<8; i++)
3228 crossterms =
new Crossterm[numCrossterms];
3229 int numRealCrossterms = 0;
3230 int *atomid = plgCterms;
3231 for(
int i=0; i<numCrossterms; i++, atomid+=8) {
3232 Crossterm *thisCrossterm = crossterms + numRealCrossterms;
3234 for(
int j=0; j<8; j++) {
3235 if(atomid[j] != lastAtomIds[j]) {
3236 duplicate_bond =
FALSE;
3238 lastAtomIds[j] = atomid[j];
3241 if(duplicate_bond) {
3244 numRealCrossterms++;
3246 thisCrossterm->
atom1 = atomid[0]-1;
3247 thisCrossterm->
atom2 = atomid[1]-1;
3248 thisCrossterm->
atom3 = atomid[2]-1;
3249 thisCrossterm->
atom4 = atomid[3]-1;
3250 thisCrossterm->
atom5 = atomid[4]-1;
3251 thisCrossterm->
atom6 = atomid[5]-1;
3252 thisCrossterm->
atom7 = atomid[6]-1;
3253 thisCrossterm->
atom8 = atomid[7]-1;
3256 atomNames[atomid[0]-1].atomtype,
3257 atomNames[atomid[1]-1].atomtype,
3258 atomNames[atomid[2]-1].atomtype,
3259 atomNames[atomid[3]-1].atomtype,
3260 atomNames[atomid[4]-1].atomtype,
3261 atomNames[atomid[5]-1].atomtype,
3262 atomNames[atomid[6]-1].atomtype,
3263 atomNames[atomid[7]-1].atomtype,
3267 numCrossterms = numRealCrossterms;
3271 occupancy =
new float[numAtoms];
3272 for(
int i=0; i<numAtoms; i++) {
3273 occupancy[i] = atomarray[i].occupancy;
3278 bfactor =
new float[numAtoms];
3279 for(
int i=0; i<numAtoms; i++) {
3280 bfactor[i] = atomarray[i].bfactor;
3288 return v1.size() > v2.size();
3298 vector<vector<int> > atomBonds(numAtoms);
3299 for (i = 0; i < numRealBonds; i++) {
3300 int a = bonds[i].atom1;
3301 int b = bonds[i].atom2;
3302 atomBonds[a].push_back(b);
3303 atomBonds[b].push_back(a);
3306 vector<int> atomsInMolecules(numAtoms, -1);
3308 for (i = 0; i < numAtoms; i++) {
3309 if (atomsInMolecules[i] == -1) {
3310 vector<int> atomList, neighborList;
3311 atomList.push_back(i);
3312 neighborList.push_back(0);
3313 int currentMolecule = molNumber++;
3315 while (atomList.size() > 0) {
3316 int particle = atomList.back();
3317 atomsInMolecules[particle] = currentMolecule;
3318 int& neighbor = neighborList.back();
3319 while (neighbor < atomBonds[particle].size() && atomsInMolecules[atomBonds[particle][neighbor]] != -1) {
3323 if (neighbor < atomBonds[particle].size()) {
3324 atomList.push_back(atomBonds[particle][neighbor]);
3325 neighborList.push_back(0);
3327 atomList.pop_back();
3328 neighborList.pop_back();
3334 numMolecules = molNumber;
3336 vector<vector<int> > molecules(numMolecules);
3337 for (i = 0; i < numAtoms; i++) {
3338 molecules[atomsInMolecules[i]].push_back(i);
3343 sort(molecules.begin(), molecules.end(),
sizeColumn);
3346 numLargeMolecules = 0;
3347 for (i = 0; i < numMolecules; i++){
3349 ++numLargeMolecules;
3361 moleculeStartIndex =
new int32[numMolecules + 1];
3362 moleculeAtom =
new int32[numAtoms];
3364 for (i = 0; i < numMolecules; i++) {
3365 moleculeStartIndex[i] = index;
3366 for (j = 0; j < molecules[i].size(); j++) {
3367 moleculeAtom[index++] = molecules[i][j];
3371 moleculeStartIndex[numMolecules] = index;
3374 printf(
"Number of Molecule with atom greater than %5d: %5d \n: ",
LARGEMOLTH, numLargeMolecules);
3375 for (i = 0; i < numMolecules; i++) {
3376 int startIdx = moleculeStartIndex[i];
3377 int endIdx = moleculeStartIndex[i+1];
3378 printf(
"Molecule Idx %5d: ", i);
3379 for (j = startIdx; j < endIdx; j++) {
3380 int atom = moleculeAtom[j];
3381 printf(
"%7d ", atom);
3399 void Molecule::build_lists_by_atom()
3403 register int numFixedAtoms = this->numFixedAtoms;
3406 if (
simParams->fixedAtomsForces ) numFixedAtoms = 0;
3413 bondsWithAtom =
new int32 *[numAtoms];
3414 cluster =
new int32 [numAtoms];
3415 clusterSize =
new int32 [numAtoms];
3417 bondsByAtom =
new int32 *[numAtoms];
3418 anglesByAtom =
new int32 *[numAtoms];
3419 dihedralsByAtom =
new int32 *[numAtoms];
3420 impropersByAtom =
new int32 *[numAtoms];
3421 crosstermsByAtom =
new int32 *[numAtoms];
3423 exclusionsByAtom =
new int32 *[numAtoms];
3424 fullExclusionsByAtom =
new int32 *[numAtoms];
3425 modExclusionsByAtom =
new int32 *[numAtoms];
3428 gromacsPairByAtom =
new int32 *[numAtoms];
3433 const int pair_self =
3436 DebugM(3,
"Building bond lists.\n");
3439 for (i=0; i<numAtoms; i++)
3443 for (i=0; i<numRealBonds; i++)
3445 byAtomSize[bonds[i].atom1]++;
3446 byAtomSize[bonds[i].atom2]++;
3448 for (i=0; i<numAtoms; i++)
3450 bondsWithAtom[i] = tmpArena->getNewArray(byAtomSize[i]+1);
3451 bondsWithAtom[i][byAtomSize[i]] = -1;
3454 for (i=0; i<numRealBonds; i++)
3456 int a1 = bonds[i].atom1;
3457 int a2 = bonds[i].atom2;
3458 bondsWithAtom[a1][byAtomSize[a1]++] = i;
3459 bondsWithAtom[a2][byAtomSize[a2]++] = i;
3467 DebugM(3,
"Calculating exclusions for QM simulation.\n");
3470 delete_qm_bonded() ;
3472 DebugM(3,
"Re-Building bond lists.\n");
3476 for (i=0; i<numAtoms; i++)
3480 for (i=0; i<numRealBonds; i++)
3482 byAtomSize[bonds[i].atom1]++;
3483 byAtomSize[bonds[i].atom2]++;
3485 for (i=0; i<numAtoms; i++)
3487 bondsWithAtom[i][byAtomSize[i]] = -1;
3490 for (i=0; i<numRealBonds; i++)
3492 int a1 = bonds[i].atom1;
3493 int a2 = bonds[i].atom2;
3494 bondsWithAtom[a1][byAtomSize[a1]++] = i;
3495 bondsWithAtom[a2][byAtomSize[a2]++] = i;
3500 for (i=0; i<numAtoms; i++) {
3503 for (i=0; i<numAtoms; i++) {
3505 while ( cluster[ci] != ci ) ci = cluster[ci];
3506 for (
int32 *b = bondsWithAtom[i]; *b != -1; ++b ) {
3507 int a = bonds[*b].atom1;
3508 if ( a == i ) a = bonds[*b].atom2;
3511 while ( cluster[ca] != ca ) ca = cluster[ca];
3512 if ( ca > ci ) cluster[ca] = cluster[ci];
3513 else cluster[ci] = cluster[ca];
3519 for (i=0; i<numAtoms; i++) {
3520 int ci = cluster[i];
3521 if ( cluster[ci] != ci ) {
3523 cluster[i] = cluster[ci];
3529 for (i=0; i<numAtoms; i++) {
3532 for (i=0; i<numAtoms; i++) {
3533 clusterSize[cluster[i]] += 1;
3546 for (i=0; i<numAtoms; i++)
3551 for (i=0; i<numBonds; i++)
3553 if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3554 && fixedAtomFlags[bonds[i].atom2] )
continue;
3556 if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1)
continue;
3557 byAtomSize[bonds[i].atom1]++;
3560 for (i=0; i<numAtoms; i++)
3562 bondsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3563 bondsByAtom[i][byAtomSize[i]] = -1;
3566 for (i=0; i<numBonds; i++)
3568 if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3569 && fixedAtomFlags[bonds[i].atom2] )
continue;
3570 if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1)
continue;
3571 int a1 = bonds[i].atom1;
3572 bondsByAtom[a1][byAtomSize[a1]++] = i;
3574 for (i=0; i<numBonds; ++i) {
3575 int a1 = bonds[i].atom1;
3576 int a2 = bonds[i].atom2;
3580 sprintf(buff,
"Atom %d is bonded to itself", a1+1);
3583 for ( j = 0; j < byAtomSize[a1]; ++j ) {
3584 int b = bondsByAtom[a1][j];
3585 int ba1 = bonds[b].atom1;
3586 int ba2 = bonds[b].atom2;
3587 if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3589 sprintf(buff,
"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3593 for ( j = 0; j < byAtomSize[a2]; ++j ) {
3594 int b = bondsByAtom[a2][j];
3595 int ba1 = bonds[b].atom1;
3596 int ba2 = bonds[b].atom2;
3597 if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3599 sprintf(buff,
"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3605 DebugM(3,
"Building angle lists.\n");
3608 for (i=0; i<numAtoms; i++)
3613 for (i=0; i<numAngles; i++)
3615 if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3616 && fixedAtomFlags[angles[i].atom2]
3617 && fixedAtomFlags[angles[i].atom3] )
continue;
3618 if ( pair_self && fepAtomFlags[angles[i].atom1] != 1)
continue;
3619 byAtomSize[angles[i].atom1]++;
3622 for (i=0; i<numAtoms; i++)
3624 anglesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3625 anglesByAtom[i][byAtomSize[i]] = -1;
3628 for (i=0; i<numAngles; i++)
3630 if ( pair_self && fepAtomFlags[angles[i].atom1] != 1)
continue;
3631 if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3632 && fixedAtomFlags[angles[i].atom2]
3633 && fixedAtomFlags[angles[i].atom3] )
continue;
3634 int a1 = angles[i].atom1;
3635 anglesByAtom[a1][byAtomSize[a1]++] = i;
3638 DebugM(3,
"Building improper lists.\n");
3641 for (i=0; i<numAtoms; i++)
3645 numCalcImpropers = 0;
3646 for (i=0; i<numImpropers; i++)
3648 if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3649 && fixedAtomFlags[impropers[i].atom2]
3650 && fixedAtomFlags[impropers[i].atom3]
3651 && fixedAtomFlags[impropers[i].atom4] )
continue;
3652 if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1)
continue;
3653 byAtomSize[impropers[i].atom1]++;
3656 for (i=0; i<numAtoms; i++)
3658 impropersByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3659 impropersByAtom[i][byAtomSize[i]] = -1;
3662 for (i=0; i<numImpropers; i++)
3664 if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3665 && fixedAtomFlags[impropers[i].atom2]
3666 && fixedAtomFlags[impropers[i].atom3]
3667 && fixedAtomFlags[impropers[i].atom4] )
continue;
3668 if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1)
continue;
3669 int a1 = impropers[i].atom1;
3670 impropersByAtom[a1][byAtomSize[a1]++] = i;
3673 DebugM(3,
"Building dihedral lists.\n");
3676 for (i=0; i<numAtoms; i++)
3680 numCalcDihedrals = 0;
3681 for (i=0; i<numDihedrals; i++)
3683 if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3684 && fixedAtomFlags[dihedrals[i].atom2]
3685 && fixedAtomFlags[dihedrals[i].atom3]
3686 && fixedAtomFlags[dihedrals[i].atom4] )
continue;
3687 if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1)
continue;
3688 byAtomSize[dihedrals[i].atom1]++;
3691 for (i=0; i<numAtoms; i++)
3693 dihedralsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3694 dihedralsByAtom[i][byAtomSize[i]] = -1;
3697 for (i=0; i<numDihedrals; i++)
3699 if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3700 && fixedAtomFlags[dihedrals[i].atom2]
3701 && fixedAtomFlags[dihedrals[i].atom3]
3702 && fixedAtomFlags[dihedrals[i].atom4] )
continue;
3703 if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1)
continue;
3704 int a1 = dihedrals[i].atom1;
3705 dihedralsByAtom[a1][byAtomSize[a1]++] = i;
3708 DebugM(3,
"Building crossterm lists.\n");
3711 for (i=0; i<numAtoms; i++)
3715 numCalcCrossterms = 0;
3716 for (i=0; i<numCrossterms; i++)
3718 if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3719 && fixedAtomFlags[crossterms[i].atom2]
3720 && fixedAtomFlags[crossterms[i].atom3]
3721 && fixedAtomFlags[crossterms[i].atom4]
3722 && fixedAtomFlags[crossterms[i].atom5]
3723 && fixedAtomFlags[crossterms[i].atom6]
3724 && fixedAtomFlags[crossterms[i].atom7]
3725 && fixedAtomFlags[crossterms[i].atom8] )
continue;
3726 if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1)
continue;
3727 byAtomSize[crossterms[i].atom1]++;
3728 numCalcCrossterms++;
3730 for (i=0; i<numAtoms; i++)
3732 crosstermsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3733 crosstermsByAtom[i][byAtomSize[i]] = -1;
3736 for (i=0; i<numCrossterms; i++)
3738 if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3739 && fixedAtomFlags[crossterms[i].atom2]
3740 && fixedAtomFlags[crossterms[i].atom3]
3741 && fixedAtomFlags[crossterms[i].atom4]
3742 && fixedAtomFlags[crossterms[i].atom5]
3743 && fixedAtomFlags[crossterms[i].atom6]
3744 && fixedAtomFlags[crossterms[i].atom7]
3745 && fixedAtomFlags[crossterms[i].atom8] )
continue;
3746 if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1)
continue;
3747 int a1 = crossterms[i].atom1;
3748 crosstermsByAtom[a1][byAtomSize[a1]++] = i;
3752 if (is_lonepairs_psf) {
3754 DebugM(3,
"Initializing lone pair host index array.\n");
3755 lphostIndexes =
new int32[numAtoms];
3756 for (i = 0; i < numAtoms; i++) {
3757 lphostIndexes[i] = -1;
3759 for (i = 0; i < numLphosts; i++) {
3760 int32 index = lphosts[i].atom1;
3761 lphostIndexes[index] = i;
3767 DebugM(3,
"Building gromacsPair lists.\n");
3770 for (i=0; i<numAtoms; i++)
3777 if ( numFixedAtoms && fixedAtomFlags[
gromacsPair[i].atom1]
3778 && fixedAtomFlags[
gromacsPair[i].atom2] )
continue;
3779 if ( pair_self && fepAtomFlags[
gromacsPair[i].atom1] != 1)
continue;
3783 for (i=0; i<numAtoms; i++)
3785 gromacsPairByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3786 gromacsPairByAtom[i][byAtomSize[i]] = -1;
3791 if ( numFixedAtoms && fixedAtomFlags[
gromacsPair[i].atom1]
3792 && fixedAtomFlags[
gromacsPair[i].atom2] )
continue;
3793 if ( pair_self && fepAtomFlags[
gromacsPair[i].atom1] != 1)
continue;
3795 gromacsPairByAtom[a1][byAtomSize[a1]++] = i;
3800 DebugM(3,
"Building exclusion data.\n");
3807 delete [] bondsWithAtom; bondsWithAtom = 0;
3808 delete tmpArena; tmpArena = 0;
3810 if (exclusions != NULL)
3811 delete [] exclusions;
3814 numTotalExclusions = exclusionSet.size();
3816 iout <<
iINFO <<
"ADDED " << (numTotalExclusions - numExclusions) <<
" IMPLICIT EXCLUSIONS\n" <<
endi;
3818 exclusions =
new Exclusion[numTotalExclusions];
3820 for ( exclIter=exclIter.begin(),i=0; exclIter != exclIter.end(); exclIter++,i++ )
3822 exclusions[i] = *exclIter;
3826 exclusionSet.clear();
3828 DebugM(3,
"Building exclusion lists.\n");
3830 for (i=0; i<numAtoms; i++)
3834 numCalcExclusions = 0;
3835 numCalcFullExclusions = 0;
3836 for (i=0; i<numTotalExclusions; i++)
3838 if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3839 && fixedAtomFlags[exclusions[i].atom2] )
continue;
3840 byAtomSize[exclusions[i].atom1]++;
3841 numCalcExclusions++;
3842 if ( ! exclusions[i].modified ) numCalcFullExclusions++;
3845 for (i=0; i<numAtoms; i++)
3847 exclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3848 exclusionsByAtom[i][byAtomSize[i]] = -1;
3851 for (i=0; i<numTotalExclusions; i++)
3853 if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3854 && fixedAtomFlags[exclusions[i].atom2] )
continue;
3855 int a1 = exclusions[i].atom1;
3856 exclusionsByAtom[a1][byAtomSize[a1]++] = i;
3861 for (i=0; i<numAtoms; i++)
3867 for (i=0; i<numTotalExclusions; i++)
3869 if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3870 && fixedAtomFlags[exclusions[i].atom2] )
continue;
3871 if ( exclusions[i].modified ) {
3872 byAtomSize2[exclusions[i].atom1]++;
3873 byAtomSize2[exclusions[i].atom2]++;
3875 byAtomSize[exclusions[i].atom1]++;
3876 byAtomSize[exclusions[i].atom2]++;
3880 for (i=0; i<numAtoms; i++)
3882 fullExclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3883 fullExclusionsByAtom[i][0] = 0;
3884 modExclusionsByAtom[i] = arena->getNewArray(byAtomSize2[i]+1);
3885 modExclusionsByAtom[i][0] = 0;
3888 for (i=0; i<numTotalExclusions; i++)
3890 int a1 = exclusions[i].atom1;
3891 int a2 = exclusions[i].atom2;
3892 if ( numFixedAtoms && fixedAtomFlags[a1]
3893 && fixedAtomFlags[a2] )
continue;
3895 if ( exclusions[i].modified ) {
3896 l1 = modExclusionsByAtom[a1];
3897 l2 = modExclusionsByAtom[a2];
3899 l1 = fullExclusionsByAtom[a1];
3900 l2 = fullExclusionsByAtom[a2];
3906 if ( ! CkMyPe() &&
simParams->printExclusions ) {
3907 for (i=0; i<numAtoms; i++) {
3908 int32 *lf = fullExclusionsByAtom[i];
3909 iout <<
"EXCL " << i <<
" FULL";
3911 for (
int j = 0; j < nf; ++j ) {
3912 iout <<
" " << *(lf++);
3915 int32 *lm = modExclusionsByAtom[i];
3916 iout <<
"EXCL " << i <<
" MOD";
3918 for (
int j = 0; j < nm; ++j ) {
3919 iout <<
" " << *(lm++);
3926 if (is_drude_psf ||
simParams->drudeOn) {
3932 if (tholes != NULL)
delete[] tholes;
3936 for (i = 0; i < numTotalExclusions; i++) {
3938 if (exclusions[i].modified)
continue;
3939 int a1 = exclusions[i].atom1;
3940 int a2 = exclusions[i].atom2;
3941 if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3947 if (numTholes != 0) tholes =
new Thole[numTholes];
3954 for (i = 0; i < numTotalExclusions; i++) {
3956 if (exclusions[i].modified)
continue;
3957 int a1 = exclusions[i].atom1;
3958 int a2 = exclusions[i].atom2;
3960 if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3961 Real thsum = drudeConsts[a1].thole + drudeConsts[a2].thole;
3962 Real aprod = drudeConsts[a1].alpha * drudeConsts[a2].alpha;
3964 Real apower = (aprod <= 0 ? 0 : powf(aprod, -1.f/6));
3965 tholes[nt].atom1 = a1;
3966 tholes[nt].atom2 = a1+1;
3967 tholes[nt].atom3 = a2;
3968 tholes[nt].atom4 = a2+1;
3969 tholes[nt].aa = apower * thsum;
3970 tholes[nt].qq = c * atoms[a1+1].charge * atoms[a2+1].charge;
3976 DebugM(3,
"Building Thole correction term lists.\n");
3977 tholesByAtom =
new int32 *[numAtoms];
3979 for (i = 0; i < numAtoms; i++) {
3983 for (i = 0; i < numTholes; i++) {
3984 if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
3985 && fixedAtomFlags[tholes[i].atom2]
3986 && fixedAtomFlags[tholes[i].atom3]
3987 && fixedAtomFlags[tholes[i].atom4] )
continue;
3988 if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1)
continue;
3989 byAtomSize[tholes[i].atom1]++;
3992 for (i = 0; i < numAtoms; i++) {
3993 tholesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3994 tholesByAtom[i][byAtomSize[i]] = -1;
3997 for (i = 0; i < numTholes; i++) {
3998 if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
3999 && fixedAtomFlags[tholes[i].atom2]
4000 && fixedAtomFlags[tholes[i].atom3]
4001 && fixedAtomFlags[tholes[i].atom4] )
continue;
4002 if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1)
continue;
4003 int a1 = tholes[i].atom1;
4004 tholesByAtom[a1][byAtomSize[a1]++] = i;
4008 DebugM(3,
"Building anisotropic term lists.\n");
4009 anisosByAtom =
new int32 *[numAtoms];
4011 for (i = 0; i < numAtoms; i++) {
4015 for (i = 0; i < numAnisos; i++) {
4016 if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
4017 && fixedAtomFlags[anisos[i].atom2]
4018 && fixedAtomFlags[anisos[i].atom3]
4019 && fixedAtomFlags[anisos[i].atom4] )
continue;
4020 if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1)
continue;
4021 byAtomSize[anisos[i].atom1]++;
4024 for (i = 0; i < numAtoms; i++) {
4025 anisosByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
4026 anisosByAtom[i][byAtomSize[i]] = -1;
4029 for (i = 0; i < numAnisos; i++) {
4030 if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
4031 && fixedAtomFlags[anisos[i].atom2]
4032 && fixedAtomFlags[anisos[i].atom3]
4033 && fixedAtomFlags[anisos[i].atom4] )
continue;
4034 if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1)
continue;
4035 int a1 = anisos[i].atom1;
4036 anisosByAtom[a1][byAtomSize[a1]++] = i;
4042 delete [] byAtomSize; byAtomSize = 0;
4043 delete [] byAtomSize2; byAtomSize2 = 0;
4049 for (i=0; i<numAtoms; i++)
4051 all_exclusions[i].
min = numAtoms;
4052 all_exclusions[i].max = -1;
4054 for (i=0; i<numTotalExclusions; i++)
4057 int a1 = exclusions[i].atom1;
4058 int a2 = exclusions[i].atom2;
4059 if ( numFixedAtoms && fixedAtomFlags[a1]
4060 && fixedAtomFlags[a2] )
continue;
4061 if ( all_exclusions[a1].min > a2 ) all_exclusions[a1].min = a2;
4062 if ( all_exclusions[a2].min > a1 ) all_exclusions[a2].min = a1;
4063 if ( a2 > all_exclusions[a1].max ) all_exclusions[a1].max = a2;
4064 if ( a1 > all_exclusions[a2].max ) all_exclusions[a2].max = a1;
4068 int maxDenseAllFull = 0;
4069 int numDenseAllFull = 0;
4070 for (i=0; i<numAtoms; i++) {
4071 int iInMiddle = ( i < all_exclusions[i].max &&
4072 i > all_exclusions[i].min ) ? 1 : 0;
4073 int s = all_exclusions[i].max - all_exclusions[i].min + 1;
4074 if ( s == fullExclusionsByAtom[i][0] + iInMiddle ) {
4075 if ( s > maxDenseAllFull ) maxDenseAllFull = s;
4076 all_exclusions[i].flags = (
char*)-1;
4078 all_exclusions[i].flags = 0;
4081 char *denseFullArray = exclArena->getNewArray(maxDenseAllFull);
4082 for ( i=0; i<maxDenseAllFull; ++i ) denseFullArray[i] =
EXCHCK_FULL;
4084 int exclmem = maxDenseAllFull;
4085 int maxExclusionFlags =
simParams->maxExclusionFlags;
4086 for (i=0; i<numAtoms; i++) {
4087 int s = all_exclusions[i].max - all_exclusions[i].min + 1;
4088 if ( all_exclusions[i].max != -1 ) {
4089 if ( all_exclusions[i].flags ) {
4090 all_exclusions[i].flags = denseFullArray;
4092 }
else if ( s < maxExclusionFlags ) {
4093 char *f = all_exclusions[i].flags = exclArena->getNewArray(s);
4094 for (
int k=0; k<s; ++k ) f[k] = 0;
4097 all_exclusions[i].flags = 0;
4100 all_exclusions[i].flags = (
char*)-1;
4104 iout <<
iINFO << numTotalExclusions <<
" exclusions consume " 4105 << exclmem <<
" bytes.\n" <<
endi;
4107 <<
" atoms sharing one array.\n" <<
endi;
4109 for (i=0; i<numTotalExclusions; i++)
4111 int a1 = exclusions[i].atom1;
4112 int a2 = exclusions[i].atom2;
4113 if ( numFixedAtoms && fixedAtomFlags[a1]
4114 && fixedAtomFlags[a2] )
continue;
4115 if ( exclusions[i].modified ) {
4116 if ( all_exclusions[a1].flags )
4117 all_exclusions[a1].flags[a2-all_exclusions[a1].min] =
EXCHCK_MOD;
4118 if ( all_exclusions[a2].flags )
4119 all_exclusions[a2].flags[a1-all_exclusions[a2].min] =
EXCHCK_MOD;
4121 if ( all_exclusions[a1].flags )
4122 all_exclusions[a1].flags[a2-all_exclusions[a1].min] =
EXCHCK_FULL;
4123 if ( all_exclusions[a2].flags )
4124 all_exclusions[a2].flags[a1-all_exclusions[a2].min] =
EXCHCK_FULL;
4153 void Molecule::build_exclusions()
4161 for (i=0; i<numExclusions; i++)
4163 exclusionSet.add(exclusions[i]);
4171 switch (exclude_flag)
4198 if (is_lonepairs_psf || is_drude_psf) {
4199 build_inherited_excl(
SCALED14 == exclude_flag);
4218 void Molecule::build_inherited_excl(
int modified) {
4220 int32 *bond1, *bond2, *bond3, *bond4, *bond5;
4221 int32 i, j, mid1, mid2, mid3, mid4;
4225 for (i = 0; i < numAtoms; i++) {
4227 if (!is_drude(i) && !is_lp(i))
continue;
4231 bond1 = bondsWithAtom[i];
4235 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4236 sprintf(err_msg,
"FOUND ISOLATED %s PARTICLE %d", idescrip, i+1);
4239 if (-1 != *(bond1+1)) {
4241 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4242 sprintf(err_msg,
"FOUND MULTIPLY LINKED %s PARTICLE %d",
4248 mid1 = bonds[*bond1].atom1;
4249 if (mid1 == i) mid1 = bonds[*bond1].atom2;
4252 if (is_drude(mid1) || is_lp(mid1) || is_hydrogen(mid1)) {
4254 const char *idescrip = (is_drude(i) ?
"DRUDE" :
"LONE PAIR");
4255 sprintf(err_msg,
"PARENT ATOM %d of %s PARTICLE %d " 4256 "IS NOT HEAVY ATOM", mid1+1, idescrip, i+1);
4260 if (exclude_flag ==
NONE) {
4270 bond2 = bondsWithAtom[mid1];
4271 while (*bond2 != -1) {
4272 j = bonds[*bond2].atom1;
4273 if ((is_drude(j) || is_lp(j)) && j != mid1) {
4274 if (i < j) exclusionSet.add(
Exclusion(i, j));
4275 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4277 j = bonds[*bond2].atom2;
4278 if ((is_drude(j) || is_lp(j)) && j != mid1) {
4279 if (i < j) exclusionSet.add(
Exclusion(i, j));
4280 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4288 bond2 = bondsWithAtom[mid1];
4291 while (*bond2 != -1) {
4292 if (bonds[*bond2].atom1 == mid1) {
4293 mid2 = bonds[*bond2].atom2;
4296 mid2 = bonds[*bond2].atom1;
4306 if (exclude_flag ==
ONETWO) {
4316 bond3 = bondsWithAtom[mid2];
4317 while (*bond3 != -1) {
4318 j = bonds[*bond3].atom1;
4319 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4320 if (i < j) exclusionSet.add(
Exclusion(i, j));
4321 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4323 j = bonds[*bond3].atom2;
4324 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4325 if (i < j) exclusionSet.add(
Exclusion(i, j));
4326 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4334 bond3 = bondsWithAtom[mid2];
4337 while (*bond3 != -1) {
4339 if (bonds[*bond3].atom1 == mid2) {
4340 mid3 = bonds[*bond3].atom2;
4343 mid3 = bonds[*bond3].atom1;
4357 else if (mid3 < i) {
4363 bond4 = bondsWithAtom[mid3];
4364 while (*bond4 != -1) {
4365 j = bonds[*bond4].atom1;
4366 if ((is_drude(j) || is_lp(j)) && j != mid3) {
4367 if (i < j) exclusionSet.add(
Exclusion(i, j));
4368 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4370 j = bonds[*bond4].atom2;
4371 if ((is_drude(j) || is_lp(j)) && j != mid3) {
4372 if (i < j) exclusionSet.add(
Exclusion(i, j));
4373 else if (j < i) exclusionSet.add(
Exclusion(j, i));
4381 bond4 = bondsWithAtom[mid3];
4384 while (*bond4 != -1) {
4386 if (bonds[*bond4].atom1 == mid3) {
4387 mid4 = bonds[*bond4].atom2;
4390 mid4 = bonds[*bond4].atom1;
4400 if (is_drude(mid4) || is_lp(mid4)) {
4405 else if (mid4 < i) {
4416 int modi = modified;
4418 int amin = (mid1 < mid4 ? mid1 : mid4);
4419 int amax = (mid1 >= mid4 ? mid1 : mid4);
4431 exclusionSet.add(
Exclusion(i, mid4, modi));
4433 else if (mid4 < i) {
4434 exclusionSet.add(
Exclusion(mid4, i, modi));
4439 bond5 = bondsWithAtom[mid4];
4440 while (*bond5 != -1) {
4441 j = bonds[*bond5].atom1;
4442 if ((is_drude(j) || is_lp(j)) && j != mid4) {
4443 if (i<j) exclusionSet.add(
Exclusion(i,j,modi));
4444 else if (j<i) exclusionSet.add(
Exclusion(j,i,modi));
4446 j = bonds[*bond5].atom2;
4447 if ((is_drude(j) || is_lp(j)) && j != mid4) {
4448 if (i<j) exclusionSet.add(
Exclusion(i,j,modi));
4449 else if (j<i) exclusionSet.add(
Exclusion(j,i,modi));
4479 void Molecule::build12excl(
void)
4485 for (i=0; i<numAtoms; i++)
4487 current_val = bondsWithAtom[i];
4490 while (*current_val != -1)
4492 if (bonds[*current_val].atom1 == i)
4494 if (i<bonds[*current_val].atom2)
4496 exclusionSet.add(
Exclusion(i,bonds[*current_val].atom2));
4501 if (i<bonds[*current_val].atom1)
4503 exclusionSet.add(
Exclusion(i,bonds[*current_val].atom1));
4519 void Molecule::build13excl(
void)
4521 int32 *bond1, *bond2;
4527 for (i=0; i<numAtoms; i++)
4529 bond1 = bondsWithAtom[i];
4532 while (*bond1 != -1)
4534 if (bonds[*bond1].atom1 == i)
4536 middle_atom=bonds[*bond1].atom2;
4540 middle_atom=bonds[*bond1].atom1;
4543 bond2 = bondsWithAtom[middle_atom];
4547 while (*bond2 != -1)
4549 if (bonds[*bond2].atom1 == middle_atom)
4551 if (i < bonds[*bond2].atom2)
4553 exclusionSet.add(
Exclusion(i,bonds[*bond2].atom2));
4558 if (i < bonds[*bond2].atom1)
4560 exclusionSet.add(
Exclusion(i,bonds[*bond2].atom1));
4580 void Molecule::build14excl(
int modified)
4582 int32 *bond1, *bond2, *bond3;
4587 for (i=0; i<numAtoms; i++)
4589 if (is_drude(i) || is_lp(i))
continue;
4592 bond1 = bondsWithAtom[i];
4594 while (*bond1 != -1)
4596 if (bonds[*bond1].atom1 == i)
4598 mid1=bonds[*bond1].atom2;
4602 mid1=bonds[*bond1].atom1;
4605 bond2 = bondsWithAtom[mid1];
4608 while (*bond2 != -1)
4610 if (bonds[*bond2].atom1 == mid1)
4612 mid2 = bonds[*bond2].atom2;
4616 mid2 = bonds[*bond2].atom1;
4628 bond3=bondsWithAtom[mid2];
4631 while (*bond3 != -1)
4633 if (bonds[*bond3].atom1 == mid2)
4635 int j = bonds[*bond3].atom2;
4641 if (i < j && !is_drude(j) && !is_lp(j))
4643 exclusionSet.add(
Exclusion(i,j,modified));
4648 int j = bonds[*bond3].atom1;
4654 if (i < j && !is_drude(j) && !is_lp(j))
4656 exclusionSet.add(
Exclusion(i,j,modified));
4678 void Molecule::stripFepExcl(
void)
4684 for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
4686 int t1 = get_fep_type(exclIter->atom1);
4687 int t2 = get_fep_type(exclIter->atom2);
4688 if (t1 && t2 && t1 !=t2 && abs(t1-t2) != 2) {
4689 fepExclusionSet.
add(*exclIter);
4692 }
else if (
simParams->pairInteractionOn ) {
4693 for ( exclIter=exclIter.begin(); exclIter != exclIter.end(); exclIter++ )
4695 int ifep_type = get_fep_type(exclIter->atom1);
4696 int jfep_type = get_fep_type(exclIter->atom2);
4699 if (ifep_type != 1 || jfep_type != 1) {
4700 fepExclusionSet.
add(*exclIter);
4705 if (!(ifep_type == 1 && jfep_type == 2) &&
4706 !(ifep_type == 2 && jfep_type == 1)) {
4707 fepExclusionSet.
add(*exclIter);
4714 for ( fepIter=fepIter.begin(); fepIter != fepIter.end(); fepIter++ )
4716 exclusionSet.del(*fepIter);
4730 if ( (psf_file =
Fopen(fname,
"r")) == NULL)
4733 sprintf(err_msg,
"UNABLE TO OPEN THE COMPRESSED .psf FILE %s", fname);
4742 NAMD_die(
"The compressed psf file format is incorrect, please re-generate!\n");
4744 float psfVer = 0.0f;
4745 sscanf(buffer,
"FORMAT VERSION: %f\n", &psfVer);
4747 NAMD_die(
"The compressed psf file format is incorrect, please re-generate!\n");
4752 NAMD_die(
"UNABLE TO FIND NSEGMENTNAMES");
4753 sscanf(buffer,
"%d", &segNamePoolSize);
4755 if(segNamePoolSize!=0)
4757 for(
int i=0; i<segNamePoolSize; i++){
4759 sscanf(buffer,
"%s", strBuf);
4760 segNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4764 for(
int i=0; i<segNamePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4769 NAMD_die(
"UNABLE TO FIND NRESIDUENAMES");
4770 sscanf(buffer,
"%d", &resNamePoolSize);
4772 if(resNamePoolSize!=0)
4774 for(
int i=0; i<resNamePoolSize; i++){
4776 sscanf(buffer,
"%s", strBuf);
4777 resNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4781 for(
int i=0; i<resNamePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4786 NAMD_die(
"UNABLE TO FIND NATOMNAMES");
4787 sscanf(buffer,
"%d", &atomNamePoolSize);
4788 if(atomNamePoolSize!=0)
4790 for(
int i=0; i<atomNamePoolSize; i++){
4792 sscanf(buffer,
"%s", strBuf);
4793 atomNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4799 NAMD_die(
"UNABLE TO FIND NATOMTYPES");
4800 sscanf(buffer,
"%d", &atomTypePoolSize);
4802 if(atomTypePoolSize!=0)
4804 for(
int i=0; i<atomTypePoolSize; i++){
4806 sscanf(buffer,
"%s", strBuf);
4807 atomTypePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
4811 for(
int i=0; i<atomTypePoolSize; i++)
NAMD_read_line(psf_file, buffer);
4816 NAMD_die(
"UNABLE TO FIND NCHARGES");
4817 sscanf(buffer,
"%d", &chargePoolSize);
4818 if(chargePoolSize!=0)
4819 atomChargePool =
new Real[chargePoolSize];
4820 for(
int i=0; i<chargePoolSize; i++){
4822 sscanf(buffer,
"%f", atomChargePool+i);
4827 NAMD_die(
"UNABLE TO FIND NMASSES");
4828 sscanf(buffer,
"%d", &massPoolSize);
4830 atomMassPool =
new Real[massPoolSize];
4831 for(
int i=0; i<massPoolSize; i++){
4833 sscanf(buffer,
"%f", atomMassPool+i);
4838 NAMD_die(
"UNABLE TO FIND ATOMSIGS");
4839 sscanf(buffer,
"%d", &atomSigPoolSize);
4842 int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
4845 for(
int i=0; i<atomSigPoolSize; i++){
4849 NAMD_die(
"UNABLE TO FIND NBONDSIGS");
4850 sscanf(buffer,
"%d", &typeCnt);
4855 for(
int j=0; j<typeCnt; j++){
4857 sscanf(buffer,
"%d | %d | %d", &tmp1, &ttype, &tisReal);
4859 oneSig.offset[0] = tmp1;
4861 if(tisReal) numRealBonds++;
4867 NAMD_die(
"UNABLE TO FIND NTHETASIGS");
4868 sscanf(buffer,
"%d", &typeCnt);
4873 for(
int j=0; j<typeCnt; j++){
4875 sscanf(buffer,
"%d %d | %d | %d", &tmp1, &tmp2, &ttype, &tisReal);
4877 oneSig.offset[0] = tmp1;
4878 oneSig.offset[1] = tmp2;
4884 NAMD_die(
"UNABLE TO FIND NPHISIGS");
4885 sscanf(buffer,
"%d", &typeCnt);
4890 for(
int j=0; j<typeCnt; j++){
4892 sscanf(buffer,
"%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
4894 oneSig.offset[0] = tmp1;
4895 oneSig.offset[1] = tmp2;
4896 oneSig.offset[2] = tmp3;
4902 NAMD_die(
"UNABLE TO FIND NIMPHISIGS");
4903 sscanf(buffer,
"%d", &typeCnt);
4908 for(
int j=0; j<typeCnt; j++){
4910 sscanf(buffer,
"%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
4912 oneSig.offset[0] = tmp1;
4913 oneSig.offset[1] = tmp2;
4914 oneSig.offset[2] = tmp3;
4920 NAMD_die(
"UNABLE TO FIND NCRTERMSIGS");
4921 sscanf(buffer,
"%d", &typeCnt);
4926 for(
int j=0; j<typeCnt; j++){
4928 sscanf(buffer,
"%d %d %d %d %d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &ttype, &tisReal);
4930 oneSig.offset[0] = tmp1;
4931 oneSig.offset[1] = tmp2;
4932 oneSig.offset[2] = tmp3;
4933 oneSig.offset[3] = tmp4;
4934 oneSig.offset[4] = tmp5;
4935 oneSig.offset[5] = tmp6;
4936 oneSig.offset[6] = tmp7;
4944 NAMD_die(
"UNABLE TO FIND NEXCLSIGS");
4946 sscanf(buffer,
"%d", &exclSigPoolSize);
4948 vector<int> fullExcls;
4949 vector<int> modExcls;
4950 for(
int i=0; i<exclSigPoolSize; i++){
4952 for(
int j=0; j<fullExclCnt; j++)
4955 for(
int j=0; j<modExclCnt; j++)
4959 exclSigPool[i].setOffsets(fullExcls, modExcls);
4968 NAMD_die(
"UNABLE TO FIND NCLUSTERS");
4970 sscanf(buffer,
"%d", &numClusters);
4975 sscanf(buffer,
"%d", &numAtoms);
4979 NAMD_die(
"UNABLE TO FIND NHYDROGENGROUP");
4980 sscanf(buffer,
"%d", &numHydrogenGroups);
4984 NAMD_die(
"UNABLE TO FIND MAXHYDROGENGROUPSIZE");
4985 sscanf(buffer,
"%d", &maxHydrogenGroupSize);
4988 NAMD_die(
"UNABLE TO FIND NMIGRATIONGROUP");
4989 sscanf(buffer,
"%d", &numMigrationGroups);
4992 NAMD_die(
"UNABLE TO FIND MAXMIGRATIONGROUPSIZE");
4993 sscanf(buffer,
"%d", &maxMigrationGroupSize);
4995 int inputRigidType = -1;
4998 NAMD_die(
"UNABLE TO FIND RIGIDBONDTYPE");
4999 sscanf(buffer,
"%d", &inputRigidType);
5002 if(
simParams->rigidBonds != inputRigidType){
5003 const char *tmpstr[]={
"RIGID_NONE",
"RIGID_ALL",
"RIGID_WATER"};
5005 sprintf(errmsg,
"RIGIDBOND TYPE MISMATCH BETWEEN INPUT (%s) AND CURRENT RUN (%s)",
5006 tmpstr[inputRigidType], tmpstr[
simParams->rigidBonds]);
5014 NAMD_die(
"UNABLE TO FIND OCCUPANCYVALID");
5015 sscanf(buffer,
"%d", &isOccupancyValid);
5018 NAMD_die(
"UNABLE TO FIND TEMPFACTORVALID");
5019 sscanf(buffer,
"%d", &isBFactorValid);
5027 build_extra_bonds(params, cfgList->
find(
"extraBondsFile"));
5031 NAMD_die(
"UNABLE TO FIND DIHEDRALPARAMARRAY");
5040 NAMD_die(
"UNABLE TO FIND IMPROPERPARAMARRAY");
5057 void Molecule::read_binary_atom_info(
int fromAtomID,
int toAtomID,
InputAtomList& inAtoms){
5058 int numAtomsPar = toAtomID-fromAtomID+1;
5059 CmiAssert(numAtomsPar > 0);
5060 CmiAssert(inAtoms.
size() == numAtomsPar);
5079 eachAtomMass = NULL;
5080 eachAtomCharge = NULL;
5082 eachAtomExclSig = NULL;
5105 FILE *perAtomFile = fopen(
simParams->binAtomFile,
"rb");
5106 if (perAtomFile==NULL) {
5108 sprintf(err_msg,
"UNABLE TO OPEN THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s",
simParams->binAtomFile);
5114 flipNum((
char *)&rMagicNum,
sizeof(
int), 1);
5116 fread(&fMagicNum,
sizeof(
int), 1, perAtomFile);
5117 if (fMagicNum==magicNum) {
5119 }
else if (fMagicNum==rMagicNum) {
5123 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS CORRUPTED",
simParams->binAtomFile);
5127 float verNum = 0.0f;
5128 fread(&verNum,
sizeof(
float), 1, perAtomFile);
5129 if (needFlip)
flipNum((
char *)&verNum,
sizeof(
float), 1);
5132 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM FILE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n",
simParams->binAtomFile);
5137 fread(&recSize,
sizeof(
int), 1, perAtomFile);
5138 if(needFlip)
flipNum((
char *)&recSize,
sizeof(
int), 1);
5141 sprintf(err_msg,
"THE ASSOCIATED PER-ATOM RECORD SIZE FOR THE COMPRESSED .psf FILE %s IS INCORRECT, PLEASE RE-GENERATE!\n",
simParams->binAtomFile);
5145 const int BUFELEMS = 32*1024;
5150 if ( _fseeki64(perAtomFile,startbyte,SEEK_CUR) )
5152 if ( fseeko(perAtomFile,startbyte,SEEK_CUR) )
5156 sprintf(errmsg,
"Error on seeking binary file %s",
simParams->binAtomFile);
5162 int atomsCnt = numAtomsPar;
5165 while(atomsCnt >= BUFELEMS) {
5166 if ( fread((
char *)elemsBuf,
sizeof(
OutputAtomRecord), BUFELEMS, perAtomFile) != BUFELEMS ) {
5168 sprintf(errmsg,
"Error on reading binary file %s",
simParams->binAtomFile);
5172 for(
int i=0; i<BUFELEMS; i++, curIdx++, oneRec++) {
5174 int aid = curIdx+fromAtomID;
5175 if(needFlip) oneRec->
flip();
5176 load_one_inputatom(aid, oneRec, fAtom);
5178 atomsCnt -= BUFELEMS;
5181 if ( fread(elemsBuf,
sizeof(
OutputAtomRecord), atomsCnt, perAtomFile) != atomsCnt ) {
5183 sprintf(errmsg,
"Error on reading binary file %s",
simParams->binAtomFile);
5187 for(
int i=curIdx; i<numAtomsPar; i++, oneRec++) {
5189 int aid = i+fromAtomID;
5190 if(needFlip) oneRec->
flip();
5191 load_one_inputatom(aid,oneRec,fAtom);
5194 if ( fclose(perAtomFile) ) {
5196 sprintf(errmsg,
"Error on closing binary file %s",
simParams->binAtomFile);
5205 is_atom_fixed(fromAtomID, &listIdx);
5206 for(
int i=listIdx; i<fixedAtomsSet->size(); i++){
5207 const AtomSet one = fixedAtomsSet->item(i);
5209 int sAtomId = one.aid1>fromAtomID ? one.aid1:fromAtomID;
5210 int eAtomId = one.aid2>toAtomID? toAtomID:one.aid2;
5211 for(
int j=sAtomId; j<=eAtomId; j++)
5212 inAtoms[j-fromAtomID].atomFixed = 1;
5219 char *thisAtomName = NULL;
5281 }
else if (thisAtomMass <= 0.05) {
5283 }
else if (thisAtomMass < 1.0) {
5285 }
else if (thisAtomMass <= 3.5) {
5287 }
else if (thisAtomName[0]==
'O' &&
5288 (thisAtomMass >= 14.0) && (thisAtomMass <= 18.0)) {
5302 void Molecule::build_excl_check_signatures(){
5304 for(
int i=0; i<exclSigPoolSize; i++){
5312 int fullMin, fullMax, modMin, modMax;
5320 if(fullMin < modMin)
5321 sigChk->
min = fullMin;
5323 sigChk->
min = modMin;
5324 if(fullMax < modMax)
5325 sigChk->
max = modMax;
5327 sigChk->
max = fullMax;
5335 iout <<
iWARN <<
"an empty exclusion signature with index " 5336 << i <<
"!\n" <<
endi;
5341 sigChk->
flags =
new char[sigChk->
max-sigChk->
min+1];
5342 memset(sigChk->
flags, 0,
sizeof(
char)*(sigChk->
max-sigChk->
min+1));
5363 void Molecule::load_atom_set(
StringList *setfile,
const char *setname,
5364 int *numAtomsInSet, AtomSetList **atomsSet)
const {
5365 if(setfile == NULL) {
5367 sprintf(errmsg,
"The text input file for %s atoms is not found!", setname);
5370 FILE *ifp = fopen(setfile->
data,
"r");
5374 sprintf(errmsg,
"ERROR IN OPENING %s ATOMS FILE: %s\n", setname, setfile->
data);
5379 int numLocalAtoms = 0;
5381 AtomSetList *localAtomsSet =
new AtomSetList();
5386 bool hasDash =
false;
5387 for(
int i=0; oneline[i] && i<128; i++){
5388 if(oneline[i]==
'-') {
5394 sscanf(oneline,
"%d-%d", &(one.aid1), &(one.aid2));
5395 if(one.aid1>one.aid2 || one.aid1<0 || one.aid2<0) {
5397 sprintf(errmsg,
"The input for %s atoms is wrong: %s\n", setname, oneline);
5400 numLocalAtoms += (one.aid2-one.aid1+1);
5402 sscanf(oneline,
"%d", &(one.aid1));
5405 sprintf(errmsg,
"The input for %s atoms is wrong: %s\n", setname, oneline);
5408 one.aid2 = one.aid1;
5411 localAtomsSet->add(one);
5415 std::sort(localAtomsSet->begin(), localAtomsSet->end());
5417 *numAtomsInSet = numLocalAtoms;
5418 *atomsSet = localAtomsSet;
5421 void Molecule::load_fixed_atoms(
StringList *fixedfile){
5422 load_atom_set(fixedfile,
"FIXED", &numFixedAtoms, &fixedAtomsSet);
5425 void Molecule::load_constrained_atoms(
StringList *constrainedfile){
5426 load_atom_set(constrainedfile,
"CONSTRAINED", &numConstraints, &constrainedAtomsSet);
5429 Bool Molecule::is_atom_in_set(AtomSetList *localAtomsSet,
int aid,
int *listIdx)
const {
5430 int idx = localAtomsSet->size();
5432 int lIdx = localAtomsSet->size()-1;
5434 while(rIdx <= lIdx){
5435 int mIdx = (rIdx+lIdx)/2;
5436 const AtomSet one = localAtomsSet->item(mIdx);
5442 }
else if(aid > one.aid1){
5446 if(listIdx) *listIdx = mIdx;
5453 if(listIdx) *listIdx = mIdx;
5459 if(listIdx) *listIdx = idx;
5477 #ifdef MEM_OPT_VERSION 5478 DebugM(3,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5487 <<
"******************************************\n" \
5488 <<
"NUM NAME TYPE RES MASS CHARGE CHARGE FEP-CHARGE" \
5489 <<
"SIGMA EPSILON SIGMA14 EPSILON14\n" \
5493 for (i=0; i<numAtoms; i++)
5495 const int vdw_type =
simParams->soluteScalingOn ?
5496 ((atoms[i].vdw_type >= LJtypecount) ?
5497 ss_vdw_type[atoms[i].vdw_type-LJtypecount] : atoms[i].vdw_type) : atoms[i].vdw_type;
5498 params->
get_vdw_params(&sigma, &epsilon, &sigma14, &epsilon14, vdw_type);
5500 DebugM(3,i+1 <<
" " << atomNames[i].atomname \
5501 <<
" " << atomNames[i].atomtype <<
" " \
5502 << atomNames[i].resname <<
" " << atoms[i].mass \
5503 <<
" " << atoms[i].charge <<
" " << sigma \
5504 <<
" " << epsilon <<
" " << sigma14 \
5505 <<
" " << epsilon14 <<
"\n" \
5523 #ifdef MEM_OPT_VERSION 5524 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5530 DebugM(2,
"BOND LIST\n" <<
"********************************\n" \
5531 <<
"ATOM1 ATOM2 TYPE1 TYPE2 k x0" \
5534 for (i=0; i<numBonds; i++)
5538 DebugM(2,bonds[i].atom1+1 <<
" " \
5539 << bonds[i].atom2+1 <<
" " \
5540 << atomNames[bonds[i].atom1].atomtype <<
" " \
5541 << atomNames[bonds[i].atom2].atomtype <<
" " << k \
5542 <<
" " << x0 <<
endi);
5560 #ifdef MEM_OPT_VERSION 5561 DebugM(2,
"WARNING: this function is not availabe in memory optimized version!\n" <<
endi);
5565 DebugM(2,
"EXPLICIT EXCLUSION LIST\n" \
5566 <<
"********************************\n" \
5570 for (i=0; i<numExclusions; i++)
5572 DebugM(2,exclusions[i].atom1+1 <<
" " \
5573 << exclusions[i].atom2+1 <<
endi);
5590 #ifdef MEM_OPT_VERSION 5596 msg->
put(massPoolSize);
5597 msg->
put(massPoolSize, atomMassPool);
5599 msg->
put(chargePoolSize);
5600 msg->
put(chargePoolSize, atomChargePool);
5603 msg->
put(atomSigPoolSize);
5604 for(
int i=0; i<atomSigPoolSize; i++)
5608 msg->
put(exclSigPoolSize);
5609 for(
int i=0; i<exclSigPoolSize; i++)
5610 exclSigPool[i].pack(msg);
5612 msg->
put(numHydrogenGroups);
5613 msg->
put(maxHydrogenGroupSize);
5614 msg->
put(numMigrationGroups);
5615 msg->
put(maxMigrationGroupSize);
5616 msg->
put(isOccupancyValid);
5617 msg->
put(isBFactorValid);
5620 msg->
put(atomNamePoolSize);
5621 for(
int i=0; i<atomNamePoolSize;i++) {
5628 int numFixedAtomsSet = fixedAtomsSet->size();
5629 msg->
put(numFixedAtoms);
5630 msg->
put(numFixedAtomsSet);
5631 msg->
put(numFixedAtomsSet*
sizeof(AtomSet), (
char *)(fixedAtomsSet->begin()));
5635 int numConstrainedAtomsSet = constrainedAtomsSet->size();
5636 msg->
put(numConstraints);
5637 msg->
put(numConstrainedAtomsSet);
5638 msg->
put(numConstrainedAtomsSet*
sizeof(AtomSet), (
char *)(constrainedAtomsSet->begin()));
5643 msg->
put(numAtoms*
sizeof(
Atom), (
char*)atoms);
5646 msg->
put(numRealBonds);
5651 msg->
put(numBonds*
sizeof(
Bond), (
char*)bonds);
5655 msg->
put(numAngles);
5658 msg->
put(numAngles*
sizeof(
Angle), (
char*)angles);
5662 msg->
put(numDihedrals);
5665 msg->
put(numDihedrals*
sizeof(
Dihedral), (
char*)dihedrals);
5669 msg->
put(num_alch_unpert_Bonds);
5670 msg->
put(num_alch_unpert_Bonds*
sizeof(
Bond), (
char*)alch_unpert_bonds);
5672 msg->
put(num_alch_unpert_Angles);
5673 msg->
put(num_alch_unpert_Angles*
sizeof(
Angle), (
char*)alch_unpert_angles);
5675 msg->
put(num_alch_unpert_Dihedrals);
5676 msg->
put(num_alch_unpert_Dihedrals*
sizeof(
Dihedral), (
char*)alch_unpert_dihedrals);
5680 msg->
put(numImpropers);
5683 msg->
put(numImpropers*
sizeof(
Improper), (
char*)impropers);
5687 msg->
put(numCrossterms);
5690 msg->
put(numCrossterms*
sizeof(
Crossterm), (
char*)crossterms);
5694 msg->
put(numDonors);
5697 msg->
put(numDonors*
sizeof(
Bond), (
char*)donors);
5701 msg->
put(numAcceptors);
5704 msg->
put(numAcceptors*
sizeof(
Bond), (
char*)acceptors);
5708 msg->
put(numExclusions);
5711 msg->
put(numExclusions*
sizeof(
Exclusion), (
char*)exclusions);
5716 msg->
put(numConstraints);
5718 msg->
put(numAtoms, consIndexes);
5722 msg->
put(numConstraints*
sizeof(ConstraintParams), (
char*)consParams);
5731 DebugM(3,
"Sending gridforce info\n" <<
endi);
5732 msg->
put(numGridforceGrids);
5734 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
5735 msg->
put(numGridforces[gridnum]);
5736 msg->
put(numAtoms, gridfrcIndexes[gridnum]);
5737 if (numGridforces[gridnum])
5739 msg->
put(numGridforces[gridnum]*
sizeof(GridforceParams), (
char*)gridfrcParams[gridnum]);
5750 msg->
put(numStirredAtoms);
5752 msg->
put(numAtoms, stirIndexes);
5754 if (numStirredAtoms)
5757 msg->
put(numStirredAtoms*
sizeof(StirParams), (
char*)stirParams);
5764 msg->
put(numMovDrag);
5765 msg->
put(numAtoms, movDragIndexes);
5768 msg->
put(numMovDrag*
sizeof(MovDragParams), (
char*)movDragParams);
5774 msg->
put(numRotDrag);
5775 msg->
put(numAtoms, rotDragIndexes);
5778 msg->
put(numRotDrag*
sizeof(RotDragParams), (
char*)rotDragParams);
5784 msg->
put(numConsTorque);
5785 msg->
put(numAtoms, consTorqueIndexes);
5788 msg->
put(numConsTorque*
sizeof(ConsTorqueParams), (
char*)consTorqueParams);
5794 { msg->
put(numConsForce);
5795 msg->
put(numAtoms, consForceIndexes);
5797 msg->
put(numConsForce*
sizeof(
Vector), (
char*)consForce);
5802 msg->
put(numMolecules);
5803 msg->
put(numLargeMolecules);
5804 msg->
put(numAtoms, moleculeAtom);
5805 msg->
put(numMolecules+1, moleculeStartIndex);
5809 msg->
put(numExPressureAtoms);
5810 msg->
put(numAtoms, exPressureAtomFlags);
5813 #ifndef MEM_OPT_VERSION 5817 msg->
put(numAtoms, langevinParams);
5823 msg->
put(numFixedAtoms);
5824 msg->
put(numAtoms, fixedAtomFlags);
5825 msg->
put(numFixedRigidBonds);
5830 msg->
put(numAtoms, qmAtomGroup);
5831 msg->
put(qmNumQMAtoms);
5832 msg->
put(qmNumQMAtoms, qmAtmChrg);
5833 msg->
put(qmNumQMAtoms, qmAtmIndx);
5835 msg->
put(qmNumBonds);
5836 msg->
put(qmMeNumBonds);
5837 msg->
put(qmMeNumBonds, qmMeMMindx);
5838 msg->
put(qmMeNumBonds, qmMeQMGrp);
5840 msg->
put(qmNumGrps);
5841 msg->
put(qmNumGrps, qmGrpID);
5842 msg->
put(qmNumGrps, qmCustPCSizes);
5843 msg->
put(qmTotCustPCs);
5844 msg->
put(qmTotCustPCs, qmCustomPCIdxs);
5850 msg->
put(numFepInitial);
5851 msg->
put(numFepFinal);
5852 msg->
put(numAtoms*
sizeof(
char), (
char*)fepAtomFlags);
5857 msg->
put(numAtoms*
sizeof(
char), (
char*)ssAtomFlags);
5858 msg->
put(ss_num_vdw_params);
5860 msg->
put(numAtoms*
sizeof(
int), (
char*)ss_index);
5863 #ifdef OPENATOM_VERSION 5866 msg->
put(numFepInitial);
5867 msg->
put(numAtoms*
sizeof(
char), (
char*)fepAtomFlags);
5869 #endif //OPENATOM_VERSION 5872 msg->
put(is_lonepairs_psf);
5873 if (is_lonepairs_psf) {
5874 msg->
put(numLphosts);
5875 msg->
put(numLphosts*
sizeof(
Lphost), (
char*)lphosts);
5877 msg->
put(is_drude_psf);
5880 msg->
put(numTholes);
5881 msg->
put(numTholes*
sizeof(
Thole), (
char*)tholes);
5882 msg->
put(numAnisos);
5883 msg->
put(numAnisos*
sizeof(
Aniso), (
char*)anisos);
5885 msg->
put(numZeroMassAtoms);
5890 msg->
put(numAtoms, (
int*)lcpoParamType);
5901 msg->
put((numAtoms),pointerToLJBeg);
5902 msg->
put((numAtoms),pointerToLJEnd);
5912 msg->
put((numAtoms),pointerToGaussBeg);
5913 msg->
put((numAtoms),pointerToGaussEnd);
5921 #ifdef MEM_OPT_VERSION 5923 build_excl_check_signatures();
5926 numBonds = numCalcBonds = 0;
5927 numAngles = numCalcAngles = 0;
5928 numDihedrals = numCalcDihedrals = 0;
5929 numImpropers = numCalcImpropers = 0;
5930 numCrossterms = numCalcCrossterms = 0;
5931 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;
5939 build_lists_by_atom();
5959 #ifdef MEM_OPT_VERSION 5963 msg->
get(massPoolSize);
5964 if(atomMassPool)
delete [] atomMassPool;
5965 atomMassPool =
new Real[massPoolSize];
5966 msg->
get(massPoolSize, atomMassPool);
5968 msg->
get(chargePoolSize);
5969 if(atomChargePool)
delete [] atomChargePool;
5970 atomChargePool =
new Real[chargePoolSize];
5971 msg->
get(chargePoolSize, atomChargePool);
5974 msg->
get(atomSigPoolSize);
5977 for(
int i=0; i<atomSigPoolSize; i++)
5981 msg->
get(exclSigPoolSize);
5982 if(exclSigPool)
delete [] exclSigPool;
5984 for(
int i=0; i<exclSigPoolSize; i++)
5985 exclSigPool[i].unpack(msg);
5987 msg->
get(numHydrogenGroups);
5988 msg->
get(maxHydrogenGroupSize);
5989 msg->
get(numMigrationGroups);
5990 msg->
get(maxMigrationGroupSize);
5991 msg->
get(isOccupancyValid);
5992 msg->
get(isBFactorValid);
5995 msg->
get(atomNamePoolSize);
5997 for(
int i=0; i<atomNamePoolSize;i++) {
6005 int numFixedAtomsSet;
6006 msg->
get(numFixedAtoms);
6007 msg->
get(numFixedAtomsSet);
6008 fixedAtomsSet =
new AtomSetList(numFixedAtomsSet);
6009 msg->
get(numFixedAtomsSet*
sizeof(AtomSet), (
char *)(fixedAtomsSet->begin()));
6013 int numConstrainedAtomsSet;
6014 msg->
get(numConstraints);
6015 msg->
get(numConstrainedAtomsSet);
6016 constrainedAtomsSet =
new AtomSetList(numConstrainedAtomsSet);
6017 msg->
get(numConstrainedAtomsSet*
sizeof(AtomSet), (
char *)(constrainedAtomsSet->begin()));
6022 atoms=
new Atom[numAtoms];
6023 msg->
get(numAtoms*
sizeof(
Atom), (
char*)atoms);
6026 msg->
get(numRealBonds);
6031 bonds=
new Bond[numBonds];
6032 msg->
get(numBonds*
sizeof(
Bond), (
char*)bonds);
6036 msg->
get(numAngles);
6040 angles=
new Angle[numAngles];
6041 msg->
get(numAngles*
sizeof(
Angle), (
char*)angles);
6045 msg->
get(numDihedrals);
6048 delete [] dihedrals;
6049 dihedrals=
new Dihedral[numDihedrals];
6050 msg->
get(numDihedrals*
sizeof(
Dihedral), (
char*)dihedrals);
6054 msg->
get(num_alch_unpert_Bonds);
6055 alch_unpert_bonds=
new Bond[num_alch_unpert_Bonds];
6056 msg->
get(num_alch_unpert_Bonds*
sizeof(
Bond), (
char*)alch_unpert_bonds);
6058 msg->
get(num_alch_unpert_Angles);
6059 alch_unpert_angles=
new Angle[num_alch_unpert_Angles];
6060 msg->
get(num_alch_unpert_Angles*
sizeof(
Angle), (
char*)alch_unpert_angles);
6062 msg->
get(num_alch_unpert_Dihedrals);
6063 alch_unpert_dihedrals=
new Dihedral[num_alch_unpert_Dihedrals];
6064 msg->
get(num_alch_unpert_Dihedrals*
sizeof(
Dihedral), (
char*)alch_unpert_dihedrals);
6068 msg->
get(numImpropers);
6071 delete [] impropers;
6072 impropers=
new Improper[numImpropers];
6073 msg->
get(numImpropers*
sizeof(
Improper), (
char*)impropers);
6077 msg->
get(numCrossterms);
6080 delete [] crossterms;
6081 crossterms=
new Crossterm[numCrossterms];
6082 msg->
get(numCrossterms*
sizeof(
Crossterm), (
char*)crossterms);
6086 msg->
get(numDonors);
6090 donors=
new Bond[numDonors];
6091 msg->
get(numDonors*
sizeof(
Bond), (
char*)donors);
6095 msg->
get(numAcceptors);
6098 delete [] acceptors;
6099 acceptors=
new Bond[numAcceptors];
6100 msg->
get(numAcceptors*
sizeof(
Bond), (
char*)acceptors);
6104 msg->
get(numExclusions);
6107 delete [] exclusions;
6108 exclusions=
new Exclusion[numExclusions];
6109 msg->
get(numExclusions*
sizeof(
Exclusion), (
char*)exclusions);
6115 msg->
get(numConstraints);
6117 delete [] consIndexes;
6118 consIndexes =
new int32[numAtoms];
6120 msg->
get(numAtoms, consIndexes);
6124 delete [] consParams;
6125 consParams =
new ConstraintParams[numConstraints];
6127 msg->
get(numConstraints*
sizeof(ConstraintParams), (
char*)consParams);
6135 DebugM(3,
"Receiving gridforce info\n");
6137 msg->
get(numGridforceGrids);
6139 DebugM(3,
"numGridforceGrids = " << numGridforceGrids <<
"\n");
6141 delete [] numGridforces;
6142 numGridforces =
new int[numGridforceGrids];
6144 delete [] gridfrcIndexes;
6145 delete [] gridfrcParams;
6146 delete [] gridfrcGrid;
6147 gridfrcIndexes =
new int32*[numGridforceGrids];
6148 gridfrcParams =
new GridforceParams*[numGridforceGrids];
6151 int grandTotalGrids = 0;
6152 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6153 msg->
get(numGridforces[gridnum]);
6155 gridfrcIndexes[gridnum] =
new int32[numAtoms];
6156 msg->
get(numAtoms, gridfrcIndexes[gridnum]);
6158 if (numGridforces[gridnum])
6160 gridfrcParams[gridnum] =
new GridforceParams[numGridforces[gridnum]];
6161 msg->
get(numGridforces[gridnum]*
sizeof(GridforceParams), (
char*)gridfrcParams[gridnum]);
6174 msg->
get(numStirredAtoms);
6176 delete [] stirIndexes;
6177 stirIndexes =
new int32[numAtoms];
6179 msg->
get(numAtoms, stirIndexes);
6181 if (numStirredAtoms)
6183 delete [] stirParams;
6184 stirParams =
new StirParams[numStirredAtoms];
6186 msg->
get(numStirredAtoms*
sizeof(StirParams), (
char*)stirParams);
6192 msg->
get(numMovDrag);
6193 delete [] movDragIndexes;
6194 movDragIndexes =
new int32[numAtoms];
6195 msg->
get(numAtoms, movDragIndexes);
6198 delete [] movDragParams;
6199 movDragParams =
new MovDragParams[numMovDrag];
6200 msg->
get(numMovDrag*
sizeof(MovDragParams), (
char*)movDragParams);
6206 msg->
get(numRotDrag);
6207 delete [] rotDragIndexes;
6208 rotDragIndexes =
new int32[numAtoms];
6209 msg->
get(numAtoms, rotDragIndexes);
6212 delete [] rotDragParams;
6213 rotDragParams =
new RotDragParams[numRotDrag];
6214 msg->
get(numRotDrag*
sizeof(RotDragParams), (
char*)rotDragParams);
6220 msg->
get(numConsTorque);
6221 delete [] consTorqueIndexes;
6222 consTorqueIndexes =
new int32[numAtoms];
6223 msg->
get(numAtoms, consTorqueIndexes);
6226 delete [] consTorqueParams;
6227 consTorqueParams =
new ConsTorqueParams[numConsTorque];
6228 msg->
get(numConsTorque*
sizeof(ConsTorqueParams), (
char*)consTorqueParams);
6234 { msg->
get(numConsForce);
6235 delete [] consForceIndexes;
6236 consForceIndexes =
new int32[numAtoms];
6237 msg->
get(numAtoms, consForceIndexes);
6239 {
delete [] consForce;
6240 consForce =
new Vector[numConsForce];
6241 msg->
get(numConsForce*
sizeof(
Vector), (
char*)consForce);
6247 msg->
get(numMolecules);
6248 msg->
get(numLargeMolecules);
6249 delete [] moleculeAtom;
6250 delete [] moleculeStartIndex;
6251 moleculeAtom =
new int32[numAtoms];
6252 moleculeStartIndex =
new int32[numMolecules+1];
6253 msg->
get(numAtoms, moleculeAtom);
6254 msg->
get(numMolecules+1, moleculeStartIndex);
6258 exPressureAtomFlags =
new int32[numAtoms];
6259 msg->
get(numExPressureAtoms);
6260 msg->
get(numAtoms, exPressureAtomFlags);
6263 #ifndef MEM_OPT_VERSION 6267 delete [] langevinParams;
6268 langevinParams =
new Real[numAtoms];
6270 msg->
get(numAtoms, langevinParams);
6276 delete [] fixedAtomFlags;
6277 fixedAtomFlags =
new int32[numAtoms];
6279 msg->
get(numFixedAtoms);
6280 msg->
get(numAtoms, fixedAtomFlags);
6281 msg->
get(numFixedRigidBonds);
6286 if( qmAtomGroup != 0)
6287 delete [] qmAtomGroup;
6288 qmAtomGroup =
new Real[numAtoms];
6290 msg->
get(numAtoms, qmAtomGroup);
6292 msg->
get(qmNumQMAtoms);
6295 delete [] qmAtmChrg;
6296 qmAtmChrg =
new Real[qmNumQMAtoms];
6298 msg->
get(qmNumQMAtoms, qmAtmChrg);
6301 delete [] qmAtmIndx;
6302 qmAtmIndx =
new int[qmNumQMAtoms];
6304 msg->
get(qmNumQMAtoms, qmAtmIndx);
6308 msg->
get(qmNumBonds);
6310 msg->
get(qmMeNumBonds);
6312 if( qmMeMMindx != 0)
6313 delete [] qmMeMMindx;
6314 qmMeMMindx =
new int[qmMeNumBonds];
6316 msg->
get(qmMeNumBonds, qmMeMMindx);
6319 delete [] qmMeQMGrp;
6320 qmMeQMGrp =
new Real[qmMeNumBonds];
6322 msg->
get(qmMeNumBonds, qmMeQMGrp);
6326 msg->
get(qmNumGrps);
6330 qmGrpID =
new Real[qmNumGrps];
6331 msg->
get(qmNumGrps, qmGrpID);
6333 if( qmCustPCSizes != 0)
6334 delete [] qmCustPCSizes;
6335 qmCustPCSizes =
new int[qmNumGrps];
6336 msg->
get(qmNumGrps, qmCustPCSizes);
6338 msg->
get(qmTotCustPCs);
6340 if( qmCustomPCIdxs != 0)
6341 delete [] qmCustomPCIdxs;
6342 qmCustomPCIdxs =
new int[qmTotCustPCs];
6343 msg->
get(qmTotCustPCs, qmCustomPCIdxs);
6349 delete [] fepAtomFlags;
6350 fepAtomFlags =
new unsigned char[numAtoms];
6352 msg->
get(numFepInitial);
6353 msg->
get(numFepFinal);
6354 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)fepAtomFlags);
6360 delete [] ssAtomFlags;
6361 delete [] ss_vdw_type;
6363 ssAtomFlags =
new unsigned char[numAtoms];
6365 ss_index =
new int [numAtoms];
6366 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)ssAtomFlags);
6367 msg->
get(ss_num_vdw_params);
6369 msg->
get(numAtoms*
sizeof(
int), (
char*)ss_index);
6372 #ifdef OPENATOM_VERSION 6375 delete [] fepAtomFlags;
6376 fepAtomFlags =
new unsigned char[numAtoms];
6378 msg->
get(numFepInitial);
6379 msg->
get(numAtoms*
sizeof(
unsigned char), (
char*)fepAtomFlags);
6380 #endif //OPENATOM_VERSION 6383 msg->
get(is_lonepairs_psf);
6384 if (is_lonepairs_psf) {
6385 msg->
get(numLphosts);
6387 lphosts =
new Lphost[numLphosts];
6388 msg->
get(numLphosts*
sizeof(
Lphost), (
char*)lphosts);
6390 msg->
get(is_drude_psf);
6392 delete[] drudeConsts;
6395 msg->
get(numTholes);
6397 tholes =
new Thole[numTholes];
6398 msg->
get(numTholes*
sizeof(
Thole), (
char*)tholes);
6399 msg->
get(numAnisos);
6401 anisos =
new Aniso[numAnisos];
6402 msg->
get(numAnisos*
sizeof(
Aniso), (
char*)anisos);
6404 msg->
get(numZeroMassAtoms);
6409 delete [] lcpoParamType;
6410 lcpoParamType =
new int[numAtoms];
6411 msg->
get(numAtoms, (
int*)lcpoParamType);
6430 delete [] gromacsPair_type;
6433 delete [] pointerToLJBeg;
6434 pointerToLJBeg =
new int[numAtoms];
6435 msg->
get((numAtoms),pointerToLJBeg);
6436 delete [] pointerToLJEnd;
6437 pointerToLJEnd =
new int[numAtoms];
6438 msg->
get((numAtoms),pointerToLJEnd);
6451 delete [] indxGaussA;
6454 delete [] indxGaussB;
6472 delete [] gRepulsive;
6475 delete [] pointerToGaussBeg;
6476 pointerToGaussBeg =
new int[numAtoms];
6477 msg->
get((numAtoms),pointerToGaussBeg);
6478 delete [] pointerToGaussEnd;
6479 pointerToGaussEnd =
new int[numAtoms];
6480 msg->
get((numAtoms),pointerToGaussEnd);
6488 #ifdef MEM_OPT_VERSION 6490 build_excl_check_signatures();
6493 numBonds = numCalcBonds = 0;
6494 numAngles = numCalcAngles = 0;
6495 numDihedrals = numCalcDihedrals = 0;
6496 numImpropers = numCalcImpropers = 0;
6497 numCrossterms = numCalcCrossterms = 0;
6498 numTotalExclusions = numCalcExclusions = numCalcFullExclusions = 0;
6506 build_atom_status();
6507 build_lists_by_atom();
6546 DebugM(3,
"Entered build_gridforce_params multi...\n");
6551 numGridforceGrids = 0;
6552 while (mgridParams != NULL) {
6553 numGridforceGrids++;
6554 mgridParams = mgridParams->
next;
6557 DebugM(3,
"numGridforceGrids = " << numGridforceGrids <<
"\n");
6558 gridfrcIndexes =
new int32*[numGridforceGrids];
6559 gridfrcParams =
new GridforceParams*[numGridforceGrids];
6561 numGridforces =
new int[numGridforceGrids];
6563 int grandTotalGrids = 0;
6565 mgridParams =
simParams->mgridforcelist.get_first();
6566 for (
int gridnum = 0; gridnum < numGridforceGrids; gridnum++) {
6567 int current_index=0;
6574 if (mgridParams == NULL) {
6575 NAMD_die(
"Problem with mgridParams!");
6581 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, gridforceFile required.");
6588 if ( (cwd == NULL) || (mgridParams->
gridforceFile[0] ==
'/') )
6594 strcpy(filename, cwd);
6598 kPDB =
new PDB(filename);
6601 NAMD_die(
"Memory allocation failed in Molecule::build_gridforce_params");
6606 NAMD_die(
"Number of atoms in grid force PDB doesn't match coordinate PDB");
6625 else if (strcasecmp(mgridParams->
gridforceCol,
"Y") == 0)
6629 else if (strcasecmp(mgridParams->
gridforceCol,
"Z") == 0)
6633 else if (strcasecmp(mgridParams->
gridforceCol,
"O") == 0)
6637 else if (strcasecmp(mgridParams->
gridforceCol,
"B") == 0)
6643 NAMD_die(
"gridforcecol must have value of X, Y, Z, O, or B");
6676 NAMD_die(
"gridforcechargecol must have value of X, Y, Z, O, or B");
6681 NAMD_die(
"gridforcecol and gridforcechargecol cannot have same value");
6688 gridfrcIndexes[gridnum] =
new int32[numAtoms];
6690 if (gridfrcIndexes[gridnum] == NULL)
6692 NAMD_die(
"memory allocation failed in Molecule::build_gridforce_params()");
6696 for (i=0; i<numAtoms; i++)
6702 kval = (kPDB->
atom(i))->xcoor();
6705 kval = (kPDB->
atom(i))->ycoor();
6708 kval = (kPDB->
atom(i))->zcoor();
6711 kval = (kPDB->
atom(i))->occupancy();
6714 kval = (kPDB->
atom(i))->temperaturefactor();
6721 gridfrcIndexes[gridnum][i] = current_index;
6727 gridfrcIndexes[gridnum][i] = -1;
6731 if (current_index == 0)
6734 iout <<
iWARN <<
"NO GRIDFORCE ATOMS WERE FOUND, BUT GRIDFORCE IS ON . . .\n" <<
endi;
6739 gridfrcParams[gridnum] =
new GridforceParams[current_index];
6740 if (gridfrcParams[gridnum] == NULL)
6742 NAMD_die(
"memory allocation failed in Molecule::build_gridforce_params");
6746 numGridforces[gridnum] = current_index;
6750 for (i=0; i<numAtoms; i++)
6752 if (gridfrcIndexes[gridnum][i] != -1)
6758 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->xcoor();
6761 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->ycoor();
6764 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->zcoor();
6767 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->occupancy();
6770 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].k = (kPDB->
atom(i))->temperaturefactor();
6778 #ifdef MEM_OPT_VERSION 6779 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atomChargePool[eachAtomCharge[i]];
6781 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = atoms[i].charge;
6785 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->xcoor();
6788 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->ycoor();
6791 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->zcoor();
6794 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->occupancy();
6797 gridfrcParams[gridnum][gridfrcIndexes[gridnum][i]].q = (kPDB->
atom(i))->temperaturefactor();
6818 strcpy(potfilename, cwd);
6825 DebugM(3,
"allocating GridforceGrid(" << gridnum <<
")\n");
6829 DebugM(4,
"grandTotalGrids = " << grandTotalGrids <<
"\n" <<
endi);
6832 mgridParams = mgridParams->
next;
6837 #ifdef MEM_OPT_VERSION 6838 void Molecule::delEachAtomSigs(){
6844 if(CmiMyRank())
return;
6846 delete [] eachAtomSig;
6847 delete [] eachAtomExclSig;
6849 eachAtomExclSig = NULL;
6852 void Molecule::delChargeSpace(){
6853 if(CmiMyRank())
return;
6855 delete [] atomChargePool;
6856 delete [] eachAtomCharge;
6857 atomChargePool = NULL;
6858 eachAtomCharge = NULL;
6861 void Molecule::delMassSpace(){
6862 if(CmiMyRank())
return;
6864 delete [] atomMassPool;
6865 delete [] eachAtomMass;
6866 atomMassPool = NULL;
6867 eachAtomMass = NULL;
6870 void Molecule::delClusterSigs() {
6871 if(CmiMyRank())
return;
6873 delete [] clusterSigs;
6877 void Molecule::delAtomNames(){
6878 if(CmiMyRank())
return;
6880 delete [] atomNames;
6885 void Molecule::delFixedAtoms(){
6886 if(CmiMyRank())
return;
6887 delete fixedAtomsSet;
6888 fixedAtomsSet = NULL;
6893 #endif // MOLECULE2_C undefined = first object file 6894 #ifdef MOLECULE2_C // second object file 6926 int current_index=0;
6936 if (consref == NULL)
6938 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, consref required.");
6939 refPDB = initial_pdb;
6943 if (consref->
next != NULL)
6945 NAMD_die(
"Multiple definitions of constraint reference file in configruation file");
6948 if ( (cwd == NULL) || (consref->
data[0] ==
'/') )
6950 strcpy(filename, consref->
data);
6954 strcpy(filename, cwd);
6955 strcat(filename, consref->
data);
6958 refPDB =
new PDB(filename);
6959 if ( refPDB == NULL )
6961 NAMD_die(
"Memory allocation failed in Molecule::build_constraint_params");
6966 NAMD_die(
"Number of atoms in constraint reference PDB doesn't match coordinate PDB");
6973 if (conskfile == NULL)
6975 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, conskfile required.");
6980 if (conskfile->
next != NULL)
6982 NAMD_die(
"Multiple definitions of constraint constant file in configuration file");
6985 if ( (consref != NULL) && (strcasecmp(consref->
data, conskfile->
data) == 0) )
6992 if ( (cwd == NULL) || (conskfile->
data[0] ==
'/') )
6994 strcpy(filename, conskfile->
data);
6998 strcpy(filename, cwd);
6999 strcat(filename, conskfile->
data);
7002 kPDB =
new PDB(filename);
7005 NAMD_die(
"Memory allocation failed in Molecule::build_constraint_params");
7010 NAMD_die(
"Number of atoms in constraint constant PDB doesn't match coordinate PDB");
7020 if (conskcol == NULL)
7026 if (conskcol->
next != NULL)
7028 NAMD_die(
"Multiple definitions of harmonic constraint column in config file");
7031 if (strcasecmp(conskcol->
data,
"X") == 0)
7035 else if (strcasecmp(conskcol->
data,
"Y") == 0)
7039 else if (strcasecmp(conskcol->
data,
"Z") == 0)
7043 else if (strcasecmp(conskcol->
data,
"O") == 0)
7047 else if (strcasecmp(conskcol->
data,
"B") == 0)
7053 NAMD_die(
"conskcol must have value of X, Y, Z, O, or B");
7060 consIndexes =
new int32[numAtoms];
7062 if (consIndexes == NULL)
7064 NAMD_die(
"memory allocation failed in Molecule::build_constraint_params()");
7068 for (i=0; i<numAtoms; i++)
7074 kval = (kPDB->
atom(i))->xcoor();
7077 kval = (kPDB->
atom(i))->ycoor();
7080 kval = (kPDB->
atom(i))->zcoor();
7083 kval = (kPDB->
atom(i))->occupancy();
7086 kval = (kPDB->
atom(i))->temperaturefactor();
7093 consIndexes[i] = current_index;
7099 consIndexes[i] = -1;
7103 if (current_index == 0)
7106 iout <<
iWARN <<
"NO CONSTRAINED ATOMS WERE FOUND, BUT CONSTRAINTS ARE ON . . .\n" <<
endi;
7111 consParams =
new ConstraintParams[current_index];
7113 if (consParams == NULL)
7115 NAMD_die(
"memory allocation failed in Molecule::build_constraint_params");
7119 numConstraints = current_index;
7123 for (i=0; i<numAtoms; i++)
7125 if (consIndexes[i] != -1)
7131 consParams[consIndexes[i]].k = (kPDB->
atom(i))->xcoor();
7134 consParams[consIndexes[i]].k = (kPDB->
atom(i))->ycoor();
7137 consParams[consIndexes[i]].k = (kPDB->
atom(i))->zcoor();
7140 consParams[consIndexes[i]].k = (kPDB->
atom(i))->occupancy();
7143 consParams[consIndexes[i]].k = (kPDB->
atom(i))->temperaturefactor();
7148 consParams[consIndexes[i]].refPos.x = (refPDB->
atom(i))->xcoor();
7149 consParams[consIndexes[i]].refPos.y = (refPDB->
atom(i))->ycoor();
7150 consParams[consIndexes[i]].refPos.z = (refPDB->
atom(i))->zcoor();
7155 if (consref != NULL)
7160 if ((conskfile != NULL) &&
7161 !((consref != NULL) &&
7162 (strcasecmp(consref->
data, conskfile->
data) == 0)
7201 int current_index=0;
7210 if (movDragFile == NULL) {
7211 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, movDragFile required.");
7216 if (movDragFile->
next != NULL) {
7217 NAMD_die(
"Multiple definitions of moving drag tag file in configuration file");
7220 if ( (cwd == NULL) || (movDragFile->
data[0] ==
'/') ) {
7221 strcpy(mainfilename, movDragFile->
data);
7223 strcpy(mainfilename, cwd);
7224 strcat(mainfilename, movDragFile->
data);
7227 tPDB =
new PDB(mainfilename);
7228 if ( tPDB == NULL ) {
7229 NAMD_die(
"Memory allocation failed in Molecule::build_movdrag_params");
7233 NAMD_die(
"Number of atoms in moving drag tag PDB doesn't match coordinate PDB");
7241 if (movDragVelFile == NULL) {
7242 if (movDragFile == NULL) {
7243 NAMD_die(
"Moving drag velocity file can not be same as coordinate PDB file");
7245 if (movDragVelFile->
next != NULL) {
7246 NAMD_die(
"Multiple definitions of moving drag velocity file in configuration file");
7253 if ( (cwd == NULL) || (movDragVelFile->
data[0] ==
'/') ) {
7254 strcpy(velfilename, movDragVelFile->
data);
7256 strcpy(velfilename, cwd);
7257 strcat(velfilename, movDragVelFile->
data);
7260 vPDB =
new PDB(velfilename);
7261 if ( vPDB == NULL ) {
7262 NAMD_die(
"Memory allocation failed in Molecule::build_movdrag_params");
7266 NAMD_die(
"Number of atoms in moving drag velocity PDB doesn't match coordinate PDB");
7278 if (movDragCol == NULL) {
7281 if (movDragCol->
next != NULL) {
7282 NAMD_die(
"Multiple definitions of drag column in config file");
7285 if (movDragFile == NULL
7286 && strcasecmp(movDragCol->
data,
"B")
7287 && strcasecmp(movDragCol->
data,
"O")) {
7288 NAMD_die(
"Can not read moving drag tags from X, Y, or Z column of the coordinate or velocity file");
7290 if (!strcasecmp(movDragCol->
data,
"X")) {
7292 }
else if (!strcasecmp(movDragCol->
data,
"Y")) {
7294 }
else if (!strcasecmp(movDragCol->
data,
"Z")) {
7296 }
else if (!strcasecmp(movDragCol->
data,
"O")) {
7298 }
else if (!strcasecmp(movDragCol->
data,
"B")) {
7302 NAMD_die(
"movDragCol must have value of X, Y, Z, O, or B");
7309 movDragIndexes =
new int32[numAtoms];
7310 if (movDragIndexes == NULL) {
7311 NAMD_die(
"memory allocation failed in Molecule::build_movdrag_params()");
7315 for (i=0; i<numAtoms; i++) {
7318 dtval = (tPDB->
atom(i))->xcoor();
7321 dtval = (tPDB->
atom(i))->ycoor();
7324 dtval = (tPDB->
atom(i))->zcoor();
7327 dtval = (tPDB->
atom(i))->occupancy();
7330 dtval = (tPDB->
atom(i))->temperaturefactor();
7336 movDragIndexes[i] = current_index;
7340 movDragIndexes[i] = -1;
7344 if (current_index == 0) {
7346 iout <<
iWARN <<
"NO DRAGGED ATOMS WERE FOUND, BUT MOVING DRAG IS ON . . . " <<
endi;
7349 movDragParams =
new MovDragParams[current_index];
7350 if (movDragParams == NULL) {
7351 NAMD_die(
"memory allocation failed in Molecule::build_movdrag_params");
7355 numMovDrag = current_index;
7359 for (i=0; i<numAtoms; i++) {
7360 if (movDragIndexes[i] != -1) {
7361 movDragParams[movDragIndexes[i]].v[0] = (vPDB->
atom(i))->xcoor();
7362 movDragParams[movDragIndexes[i]].v[1] = (vPDB->
atom(i))->ycoor();
7363 movDragParams[movDragIndexes[i]].v[2] = (vPDB->
atom(i))->zcoor();
7367 if (movDragFile != NULL)
delete tPDB;
7368 if (movDragVelFile != NULL)
delete vPDB;
7405 PDB *tPDB, *aPDB, *pPDB, *vPDB;
7407 int current_index=0;
7420 if (rotDragFile == NULL) {
7421 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, rotDragFile required.");
7426 if (rotDragFile->
next != NULL) {
7427 NAMD_die(
"Multiple definitions of rotating drag tag file in configuration file");
7430 if ( (cwd == NULL) || (rotDragFile->
data[0] ==
'/') ) {
7431 strcpy(mainfilename, rotDragFile->
data);
7433 strcpy(mainfilename, cwd);
7434 strcat(mainfilename, rotDragFile->
data);
7437 tPDB =
new PDB(mainfilename);
7438 if ( tPDB == NULL ) {
7439 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7443 NAMD_die(
"Number of atoms in rotating drag tag PDB doesn't match coordinate PDB");
7451 if (rotDragAxisFile == NULL) {
7452 if (rotDragFile == NULL) {
7453 NAMD_die(
"Rotating drag axis file can not be same as coordinate PDB file");
7455 if (rotDragAxisFile->
next != NULL) {
7456 NAMD_die(
"Multiple definitions of rotating drag axis file in configuration file");
7458 if (rotDragPivotFile == NULL) {
7459 NAMD_die(
"Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
7466 if ( (cwd == NULL) || (rotDragAxisFile->
data[0] ==
'/') ) {
7467 strcpy(axisfilename, rotDragAxisFile->
data);
7469 strcpy(axisfilename, cwd);
7470 strcat(axisfilename, rotDragAxisFile->
data);
7473 aPDB =
new PDB(axisfilename);
7474 if ( aPDB == NULL ) {
7475 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7479 NAMD_die(
"Number of atoms in rotating drag axis PDB doesn't match coordinate PDB");
7487 if (rotDragPivotFile == NULL) {
7488 if (rotDragFile == NULL) {
7489 NAMD_die(
"Rotating drag pivot point file can not be same as coordinate PDB file");
7491 if (rotDragPivotFile->
next != NULL) {
7492 NAMD_die(
"Multiple definitions of rotating drag pivot point file in configuration file");
7494 if (rotDragAxisFile == NULL) {
7495 NAMD_die(
"Need to specify at least one of rotDragAxisFile and rotDragPivotFile; they can not be same");
7502 if ( (cwd == NULL) || (rotDragPivotFile->
data[0] ==
'/') ) {
7503 strcpy(pivotfilename, rotDragPivotFile->
data);
7505 strcpy(pivotfilename, cwd);
7506 strcat(pivotfilename, rotDragPivotFile->
data);
7509 pPDB =
new PDB(pivotfilename);
7510 if ( pPDB == NULL ) {
7511 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7515 NAMD_die(
"Number of atoms in rotating drag pivot point PDB doesn't match coordinate PDB");
7524 if (rotDragVelFile == NULL) {
7527 if (rotDragVelFile->
next != NULL) {
7528 NAMD_die(
"Multiple definitions of rotating drag velocity file in configuration file");
7531 if ( (cwd == NULL) || (rotDragVelFile->
data[0] ==
'/') ) {
7532 strcpy(velfilename, rotDragVelFile->
data);
7534 strcpy(velfilename, cwd);
7535 strcat(velfilename, rotDragVelFile->
data);
7538 vPDB =
new PDB(velfilename);
7539 if ( vPDB == NULL ) {
7540 NAMD_die(
"Memory allocation failed in Molecule::build_rotdrag_params");
7544 NAMD_die(
"Number of atoms in rotating drag velocity PDB doesn't match coordinate PDB");
7555 if (rotDragCol == NULL) {
7558 if (rotDragCol->
next != NULL) {
7559 NAMD_die(
"Multiple definitions of drag tag column in config file");
7562 if ( rotDragFile == NULL
7563 && (!strcasecmp(rotDragCol->
data,
"X")
7564 || !strcasecmp(rotDragCol->
data,
"Y")
7565 || !strcasecmp(rotDragCol->
data,
"Z"))) {
7566 NAMD_die(
"Can not read rotating drag tags from X, Y, or Z column of the PDB coordinate file");
7568 if (!strcasecmp(rotDragCol->
data,
"X")) {
7570 }
else if (!strcasecmp(rotDragCol->
data,
"Y")) {
7572 }
else if (!strcasecmp(rotDragCol->
data,
"Z")) {
7574 }
else if (!strcasecmp(rotDragCol->
data,
"O")) {
7576 }
else if (!strcasecmp(rotDragCol->
data,
"B")) {
7580 NAMD_die(
"rotDragCol must have value of X, Y, Z, O, or B");
7592 if (rotDragVelCol == NULL) {
7595 if (rotDragVelCol->
next != NULL) {
7596 NAMD_die(
"Multiple definitions of drag angular velocity column in config file");
7599 if (rotDragVelFile == NULL
7600 && rotDragFile == NULL
7601 && strcasecmp(rotDragCol->
data,
"B")
7602 && strcasecmp(rotDragCol->
data,
"O")) {
7603 NAMD_die(
"Can not read rotating drag angular velocities from X, Y, or Z column of the PDB coordinate file");
7605 if (!strcasecmp(rotDragVelCol->
data,
"X")) {
7607 }
else if (!strcasecmp(rotDragVelCol->
data,
"Y")) {
7609 }
else if (!strcasecmp(rotDragVelCol->
data,
"Z")) {
7611 }
else if (!strcasecmp(rotDragVelCol->
data,
"O")) {
7613 }
else if (!strcasecmp(rotDragVelCol->
data,
"B")) {
7617 NAMD_die(
"rotDragVelCol must have value of X, Y, Z, O, or B");
7624 rotDragIndexes =
new int32[numAtoms];
7625 if (rotDragIndexes == NULL) {
7626 NAMD_die(
"memory allocation failed in Molecule::build_rotdrag_params()");
7630 for (i=0; i<numAtoms; i++) {
7633 dtval = (tPDB->
atom(i))->xcoor();
7636 dtval = (tPDB->
atom(i))->ycoor();
7639 dtval = (tPDB->
atom(i))->zcoor();
7642 dtval = (tPDB->
atom(i))->occupancy();
7645 dtval = (tPDB->
atom(i))->temperaturefactor();
7651 rotDragIndexes[i] = current_index;
7655 rotDragIndexes[i] = -1;
7659 if (current_index == 0) {
7660 iout <<
iWARN <<
"NO DRAGGED ATOMS WERE FOUND, BUT ROTATING DRAG IS ON . . . " <<
endi;
7662 rotDragParams =
new RotDragParams[current_index];
7663 if (rotDragParams == NULL) {
7664 NAMD_die(
"memory allocation failed in Molecule::build_rotdrag_params");
7668 numRotDrag = current_index;
7672 for (i=0; i<numAtoms; i++) {
7673 if (rotDragIndexes[i] != -1) {
7674 rotDragParams[rotDragIndexes[i]].a[0] = (aPDB->
atom(i))->xcoor();
7675 rotDragParams[rotDragIndexes[i]].a[1] = (aPDB->
atom(i))->ycoor();
7676 rotDragParams[rotDragIndexes[i]].a[2] = (aPDB->
atom(i))->zcoor();
7677 rotDragParams[rotDragIndexes[i]].p[0] = (pPDB->
atom(i))->xcoor();
7678 rotDragParams[rotDragIndexes[i]].p[1] = (pPDB->
atom(i))->ycoor();
7679 rotDragParams[rotDragIndexes[i]].p[2] = (pPDB->
atom(i))->zcoor();
7682 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->xcoor();
7685 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->ycoor();
7688 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->zcoor();
7691 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->occupancy();
7694 rotDragParams[rotDragIndexes[i]].v = (vPDB->
atom(i))->temperaturefactor();
7700 if (rotDragFile != NULL)
delete tPDB;
7701 if (rotDragAxisFile != NULL)
delete aPDB;
7702 if (rotDragPivotFile != NULL)
delete pPDB;
7703 if (rotDragVelFile != NULL)
delete vPDB;
7740 PDB *tPDB, *aPDB, *pPDB, *vPDB;
7742 int current_index=0;
7755 if (consTorqueFile == NULL) {
7756 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, consTorqueFile required.");
7761 if (consTorqueFile->
next != NULL) {
7762 NAMD_die(
"Multiple definitions of \"constant\" torque tag file in configuration file");
7765 if ( (cwd == NULL) || (consTorqueFile->
data[0] ==
'/') ) {
7766 strcpy(mainfilename, consTorqueFile->
data);
7768 strcpy(mainfilename, cwd);
7769 strcat(mainfilename, consTorqueFile->
data);
7772 tPDB =
new PDB(mainfilename);
7773 if ( tPDB == NULL ) {
7774 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7778 NAMD_die(
"Number of atoms in \"constant\" torque tag PDB doesn't match coordinate PDB");
7786 if (consTorqueAxisFile == NULL) {
7787 if (consTorqueFile == NULL) {
7788 NAMD_die(
"\"Constant\" torque axis file can not be same as coordinate PDB file");
7790 if (consTorqueAxisFile->
next != NULL) {
7791 NAMD_die(
"Multiple definitions of \"constant\" torque axis file in configuration file");
7793 if (consTorquePivotFile == NULL) {
7794 NAMD_die(
"Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
7801 if ( (cwd == NULL) || (consTorqueAxisFile->
data[0] ==
'/') ) {
7802 strcpy(axisfilename, consTorqueAxisFile->
data);
7804 strcpy(axisfilename, cwd);
7805 strcat(axisfilename, consTorqueAxisFile->
data);
7808 aPDB =
new PDB(axisfilename);
7809 if ( aPDB == NULL ) {
7810 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7814 NAMD_die(
"Number of atoms in \"constant\" torque axis PDB doesn't match coordinate PDB");
7822 if (consTorquePivotFile == NULL) {
7823 if (consTorqueFile == NULL) {
7824 NAMD_die(
"\"Constant\" torque pivot point file can not be same as coordinate PDB file");
7826 if (consTorquePivotFile->
next != NULL) {
7827 NAMD_die(
"Multiple definitions of \"constant\" torque pivot point file in configuration file");
7829 if (consTorqueAxisFile == NULL) {
7830 NAMD_die(
"Need to specify at least one of consTorqueAxisFile and consTorquePivotFile; they can not be same");
7837 if ( (cwd == NULL) || (consTorquePivotFile->
data[0] ==
'/') ) {
7838 strcpy(pivotfilename, consTorquePivotFile->
data);
7840 strcpy(pivotfilename, cwd);
7841 strcat(pivotfilename, consTorquePivotFile->
data);
7844 pPDB =
new PDB(pivotfilename);
7845 if ( pPDB == NULL ) {
7846 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7850 NAMD_die(
"Number of atoms in \"constant\" torque pivot point PDB doesn't match coordinate PDB");
7859 if (consTorqueValFile == NULL) {
7862 if (consTorqueValFile->
next != NULL) {
7863 NAMD_die(
"Multiple definitions of \"constant\" torque velocity file in configuration file");
7866 if ( (cwd == NULL) || (consTorqueValFile->
data[0] ==
'/') ) {
7867 strcpy(velfilename, consTorqueValFile->
data);
7869 strcpy(velfilename, cwd);
7870 strcat(velfilename, consTorqueValFile->
data);
7873 vPDB =
new PDB(velfilename);
7874 if ( vPDB == NULL ) {
7875 NAMD_die(
"Memory allocation failed in Molecule::build_constorque_params");
7879 NAMD_die(
"Number of atoms in \"constant\" torque velocity PDB doesn't match coordinate PDB");
7890 if (consTorqueCol == NULL) {
7893 if (consTorqueCol->
next != NULL) {
7894 NAMD_die(
"Multiple definitions of torque tag column in config file");
7897 if ( consTorqueFile == NULL
7898 && (!strcasecmp(consTorqueCol->
data,
"X")
7899 || !strcasecmp(consTorqueCol->
data,
"Y")
7900 || !strcasecmp(consTorqueCol->
data,
"Z"))) {
7901 NAMD_die(
"Can not read \"constant\" torque tags from X, Y, or Z column of the PDB coordinate file");
7903 if (!strcasecmp(consTorqueCol->
data,
"X")) {
7905 }
else if (!strcasecmp(consTorqueCol->
data,
"Y")) {
7907 }
else if (!strcasecmp(consTorqueCol->
data,
"Z")) {
7909 }
else if (!strcasecmp(consTorqueCol->
data,
"O")) {
7911 }
else if (!strcasecmp(consTorqueCol->
data,
"B")) {
7915 NAMD_die(
"consTorqueCol must have value of X, Y, Z, O, or B");
7927 if (consTorqueValCol == NULL) {
7930 if (consTorqueValCol->
next != NULL) {
7931 NAMD_die(
"Multiple definitions of torque value column in config file");
7934 if (consTorqueValFile == NULL
7935 && consTorqueFile == NULL
7936 && strcasecmp(consTorqueCol->
data,
"B")
7937 && strcasecmp(consTorqueCol->
data,
"O")) {
7938 NAMD_die(
"Can not read \"constant\" torque values from X, Y, or Z column of the PDB coordinate file");
7940 if (!strcasecmp(consTorqueValCol->
data,
"X")) {
7942 }
else if (!strcasecmp(consTorqueValCol->
data,
"Y")) {
7944 }
else if (!strcasecmp(consTorqueValCol->
data,
"Z")) {
7946 }
else if (!strcasecmp(consTorqueValCol->
data,
"O")) {
7948 }
else if (!strcasecmp(consTorqueValCol->
data,
"B")) {
7952 NAMD_die(
"consTorqueValCol must have value of X, Y, Z, O, or B");
7959 consTorqueIndexes =
new int32[numAtoms];
7960 if (consTorqueIndexes == NULL) {
7961 NAMD_die(
"memory allocation failed in Molecule::build_constorque_params()");
7965 for (i=0; i<numAtoms; i++) {
7968 dtval = (tPDB->
atom(i))->xcoor();
7971 dtval = (tPDB->
atom(i))->ycoor();
7974 dtval = (tPDB->
atom(i))->zcoor();
7977 dtval = (tPDB->
atom(i))->occupancy();
7980 dtval = (tPDB->
atom(i))->temperaturefactor();
7986 consTorqueIndexes[i] = current_index;
7990 consTorqueIndexes[i] = -1;
7994 if (current_index == 0) {
7995 iout <<
iWARN <<
"NO TORQUED ATOMS WERE FOUND, BUT \"CONSTANT\" TORQUE IS ON . . . " <<
endi;
7997 consTorqueParams =
new ConsTorqueParams[current_index];
7998 if (consTorqueParams == NULL) {
7999 NAMD_die(
"memory allocation failed in Molecule::build_constorque_params");
8003 numConsTorque = current_index;
8007 for (i=0; i<numAtoms; i++) {
8008 if (consTorqueIndexes[i] != -1) {
8009 consTorqueParams[consTorqueIndexes[i]].a[0] = (aPDB->
atom(i))->xcoor();
8010 consTorqueParams[consTorqueIndexes[i]].a[1] = (aPDB->
atom(i))->ycoor();
8011 consTorqueParams[consTorqueIndexes[i]].a[2] = (aPDB->
atom(i))->zcoor();
8012 consTorqueParams[consTorqueIndexes[i]].p[0] = (pPDB->
atom(i))->xcoor();
8013 consTorqueParams[consTorqueIndexes[i]].p[1] = (pPDB->
atom(i))->ycoor();
8014 consTorqueParams[consTorqueIndexes[i]].p[2] = (pPDB->
atom(i))->zcoor();
8017 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->xcoor();
8020 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->ycoor();
8023 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->zcoor();
8026 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->occupancy();
8029 consTorqueParams[consTorqueIndexes[i]].v = (vPDB->
atom(i))->temperaturefactor();
8035 if (consTorqueFile != NULL)
delete tPDB;
8036 if (consTorqueAxisFile != NULL)
delete aPDB;
8037 if (consTorquePivotFile != NULL)
delete pPDB;
8038 if (consTorqueValFile != NULL)
delete vPDB;
8064 iout <<
iWARN <<
"NO CONSTANT FORCES SPECIFIED, BUT CONSTANT FORCE IS ON . . .\n" <<
endi;
8065 consForceIndexes =
new int32[numAtoms];
8066 for (i=0; i<numAtoms; i++) consForceIndexes[i] = -1;
8070 if ((forcePDB=
new PDB(filename)) == NULL)
8071 NAMD_die(
"Memory allocation failed in Molecule::build_constant_forces");
8073 NAMD_die(
"Number of atoms in constant force PDB doesn't match coordinate PDB");
8078 consForceIndexes =
new int32[numAtoms];
8079 if (consForceIndexes == NULL)
8080 NAMD_die(
"memory allocation failed in Molecule::build_constant_forces()");
8084 for (i=0; i<numAtoms; i++)
8088 consForceIndexes[i] = -1;
8091 consForceIndexes[i] = numConsForce++;
8093 if (numConsForce == 0)
8095 iout <<
iWARN <<
"NO NON-ZERO FORCES WERE FOUND, BUT CONSTANT FORCE IS ON . . .\n" <<
endi;
8098 consForce =
new Vector[numConsForce];
8099 if (consForce == NULL)
8100 NAMD_die(
"memory allocation failed in Molecule::build_constant_forces");
8102 for (i=0; i<numAtoms; i++)
8103 if ((index=consForceIndexes[i]) != -1)
8120 langevinParams =
new Real[numAtoms];
8122 if ( (langevinParams == NULL) )
8124 NAMD_die(
"memory allocation failed in Molecule::build_langevin_params()");
8128 for (
int i=0; i<numAtoms; i++)
8132 if ( (! doHydrogen) && is_hydrogen(i) ) bval = 0;
8133 else if ( is_lp(i) ) bval = 0;
8134 else if ( is_drude(i) ) bval = drudeCoupling;
8137 langevinParams[i] = bval;
8174 if (langfile == NULL)
8176 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, langevinFile required.");
8181 if (langfile->
next != NULL)
8183 NAMD_die(
"Multiple definitions of langvein PDB file in configuration file");
8186 if ( (cwd == NULL) || (langfile->
data[0] ==
'/') )
8188 strcpy(filename, langfile->
data);
8192 strcpy(filename, cwd);
8193 strcat(filename, langfile->
data);
8196 bPDB =
new PDB(filename);
8199 NAMD_die(
"Memory allocation failed in Molecule::build_langevin_params");
8204 NAMD_die(
"Number of atoms in langevin parameter PDB doesn't match coordinate PDB");
8213 if (langcol == NULL)
8219 if (langcol->
next != NULL)
8221 NAMD_die(
"Multiple definitions of langevin parameter column in config file");
8224 if (strcasecmp(langcol->
data,
"X") == 0)
8228 else if (strcasecmp(langcol->
data,
"Y") == 0)
8232 else if (strcasecmp(langcol->
data,
"Z") == 0)
8236 else if (strcasecmp(langcol->
data,
"O") == 0)
8240 else if (strcasecmp(langcol->
data,
"B") == 0)
8246 NAMD_die(
"langevincol must have value of X, Y, Z, O, or B");
8251 langevinParams =
new Real[numAtoms];
8253 if ( (langevinParams == NULL) )
8255 NAMD_die(
"memory allocation failed in Molecule::build_langevin_params()");
8259 for (i=0; i<numAtoms; i++)
8265 bval = (bPDB->
atom(i))->xcoor();
8268 bval = (bPDB->
atom(i))->ycoor();
8271 bval = (bPDB->
atom(i))->zcoor();
8274 bval = (bPDB->
atom(i))->occupancy();
8277 bval = (bPDB->
atom(i))->temperaturefactor();
8282 langevinParams[i] = bval;
8286 if (langfile != NULL)
8325 if (fixedfile == NULL)
8327 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, fixedAtomsFile required.");
8332 if (fixedfile->
next != NULL)
8334 NAMD_die(
"Multiple definitions of fixed atoms PDB file in configuration file");
8337 if ( (cwd == NULL) || (fixedfile->
data[0] ==
'/') )
8339 strcpy(filename, fixedfile->
data);
8343 strcpy(filename, cwd);
8344 strcat(filename, fixedfile->
data);
8347 bPDB =
new PDB(filename);
8350 NAMD_die(
"Memory allocation failed in Molecule::build_fixed_atoms");
8355 NAMD_die(
"Number of atoms in fixed atoms PDB doesn't match coordinate PDB");
8364 if (fixedcol == NULL)
8370 if (fixedcol->
next != NULL)
8372 NAMD_die(
"Multiple definitions of fixed atoms column in config file");
8375 if (strcasecmp(fixedcol->
data,
"X") == 0)
8379 else if (strcasecmp(fixedcol->
data,
"Y") == 0)
8383 else if (strcasecmp(fixedcol->
data,
"Z") == 0)
8387 else if (strcasecmp(fixedcol->
data,
"O") == 0)
8391 else if (strcasecmp(fixedcol->
data,
"B") == 0)
8397 NAMD_die(
"fixedatomscol must have value of X, Y, Z, O, or B");
8402 fixedAtomFlags =
new int32[numAtoms];
8404 if (fixedAtomFlags == NULL)
8406 NAMD_die(
"memory allocation failed in Molecule::build_fixed_atoms()");
8412 for (i=0; i<numAtoms; i++)
8418 bval = (bPDB->
atom(i))->xcoor();
8421 bval = (bPDB->
atom(i))->ycoor();
8424 bval = (bPDB->
atom(i))->zcoor();
8427 bval = (bPDB->
atom(i))->occupancy();
8430 bval = (bPDB->
atom(i))->temperaturefactor();
8436 fixedAtomFlags[i] = 1;
8440 fixedAtomFlags[i] = 0;
8445 if (fixedfile != NULL)
8452 if ( numRigidBonds ) {
8454 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
8455 int parentIsFixed = 0;
8456 for( ; h_i != h_e; ++h_i ) {
8458 parentIsFixed = fixedAtomFlags[h_i->
atomID];
8459 if ( (rigidBondLengths[h_i->
atomID] > 0.)
8460 && fixedAtomFlags[h_i[1].atomID]
8461 && fixedAtomFlags[h_i[2].atomID] ) {
8462 ++numFixedRigidBonds;
8465 if ( (rigidBondLengths[h_i->
atomID] > 0.)
8466 && fixedAtomFlags[h_i->
atomID]
8467 && parentIsFixed ) {
8468 ++numFixedRigidBonds;
8478 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
8480 for( ; h_i != h_e; ++h_i ) {
8482 if ( allFixed ) ++numFixedGroups;
8485 allFixed = allFixed && fixedAtomFlags[h_i->
atomID];
8487 if ( allFixed ) ++numFixedGroups;
8524 int current_index=0;
8532 if (stirredfile == NULL)
8534 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, stirFilename required.");
8537 iout <<
iWARN <<
"STIRRING USING INITIAL POSITION FILE FOR REFERENCE POSITIONS" <<
endi;
8541 if (stirredfile->
next != NULL)
8543 NAMD_die(
"Multiple definitions of stirred atoms PDB file in configuration file");
8546 if ( (cwd == NULL) || (stirredfile->
data[0] ==
'/') )
8548 strcpy(filename, stirredfile->
data);
8552 strcpy(filename, cwd);
8553 strcat(filename, stirredfile->
data);
8557 sPDB =
new PDB(filename);
8561 NAMD_die(
"Memory allocation failed in Molecule::build_stirred_atoms");
8567 NAMD_die(
"Number of atoms in stirred atoms PDB doesn't match coordinate PDB");
8579 if (stirredcol == NULL)
8585 if (stirredcol->
next != NULL)
8587 NAMD_die(
"Multiple definitions of stirred atoms column in config file");
8590 if (strcasecmp(stirredcol->
data,
"O") == 0)
8594 else if (strcasecmp(stirredcol->
data,
"B") == 0)
8600 NAMD_die(
"stirredAtomsCol must have value of O or B");
8607 stirIndexes =
new int32[numAtoms];
8609 if (stirIndexes == NULL)
8611 NAMD_die(
"memory allocation failed in Molecule::build_stirred_params()");
8616 for (i=0; i<numAtoms; i++)
8626 bval = (sPDB->
atom(i))->occupancy();
8629 bval = (sPDB->
atom(i))->temperaturefactor();
8638 stirIndexes[i] = current_index;
8644 stirIndexes[i] = -1;
8652 if (current_index == 0)
8655 iout <<
iWARN <<
"NO STIRRED ATOMS WERE FOUND, BUT STIRRING TORQUES ARE ON . . .\n" <<
endi;
8660 stirParams =
new StirParams[current_index];
8662 if (stirParams == NULL)
8664 NAMD_die(
"memory allocation failed in Molecule::build_stir_params");
8668 numStirredAtoms = current_index;
8672 for (i=0; i<numAtoms; i++)
8674 if (stirIndexes[i] != -1)
8678 stirParams[stirIndexes[i]].refPos.x = (sPDB->
atom(i))->xcoor();
8679 stirParams[stirIndexes[i]].refPos.y = (sPDB->
atom(i))->ycoor();
8680 stirParams[stirIndexes[i]].refPos.z = (sPDB->
atom(i))->zcoor();
8685 if (stirredfile != NULL)
8701 int a1,a2,a3,a4;
float k, ref, upper;
8702 int anglesNormal = (
simParams->extraBondsCosAngles ? 0 : 1 );
8703 #ifndef MEM_OPT_VERSION 8716 NAMD_die(
"NO EXTRA BONDS FILES SPECIFIED");
8719 for ( ; file; file = file->
next ) {
8720 FILE *f = fopen(file->
data,
"r");
8722 sprintf(err_msg,
"UNABLE TO OPEN EXTRA BONDS FILE %s", file->
data);
8734 if (ret_code!=0)
break;
8737 sscanf(buffer,
"%s",type);
8739 #define CHECKATOMID(ATOMID) if ( ATOMID < 0 || ATOMID >= numAtoms ) badatom = 1; 8743 if ( ! strncasecmp(type,
"bond",4) ) {
8744 if ( sscanf(buffer,
"%s %d %d %f %f %s",
8745 type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1;
8751 #ifndef MEM_OPT_VERSION 8759 tmpv.
k = k; tmpv.
x0 = ref;
8761 }
else if ( ! strncasecmp(type,
"wall",4) ) {
8764 if ( sscanf(buffer,
"%s %d %d %f %f %f %s",
8765 type, &a1, &a2, &k, &ref, &upper, err_msg) != 6 ) badline = 1;
8766 else if (upper < ref) badline = 1;
8772 #ifndef MEM_OPT_VERSION 8780 tmpv.
k = k; tmpv.
x0 = ref; tmpv.
x1 = upper;
8782 }
else if ( ! strncasecmp(type,
"angle",4) ) {
8783 if ( sscanf(buffer,
"%s %d %d %d %f %f %s",
8784 type, &a1, &a2, &a3, &k, &ref, err_msg) != 6 ) badline = 1;
8790 #ifndef MEM_OPT_VERSION 8798 tmpv.
k = k; tmpv.
theta0 = ref / 180. *
PI;
8800 tmpv.
normal = anglesNormal;
8803 }
else if ( ! strncasecmp(type,
"dihedral",4) ) {
8805 int ret = 1 + sscanf(buffer,
"%s %d %d %d %d %f %f %s",
8806 type, &a1, &a2, &a3, &a4, &k, &ref, err_msg);
8808 ret = sscanf(buffer,
"%s %d %d %d %d %f %d %f %s",
8809 type, &a1, &a2, &a3, &a4, &k, &n, &ref, err_msg);
8811 if ( ret != 8 ) badline = 1;
8818 #ifndef MEM_OPT_VERSION 8829 }
else if ( ! strncasecmp(type,
"improper",4) ) {
8831 int ret = 1 + sscanf(buffer,
"%s %d %d %d %d %f %f %s",
8832 type, &a1, &a2, &a3, &a4, &k, &ref, err_msg);
8834 ret = sscanf(buffer,
"%s %d %d %d %d %f %d %f %s",
8835 type, &a1, &a2, &a3, &a4, &k, &n, &ref, err_msg);
8837 if ( ret != 8 ) badline = 1;
8844 #ifndef MEM_OPT_VERSION 8855 }
else if ( ! strncasecmp(type,
"#",1) ) {
8862 sprintf(err_msg,
"BAD LINE IN EXTRA BONDS FILE %s: %s",
8863 file->
data, buffer);
8867 sprintf(err_msg,
"BAD ATOM ID IN EXTRA BONDS FILE %s: %s",
8868 file->
data, buffer);
8877 if ( extraNumBonds ) {
8878 iout <<
iINFO <<
"READ " << extraNumBonds <<
" EXTRA BONDS\n" <<
endi;
8880 #ifndef MEM_OPT_VERSION 8881 Bond *newbonds =
new Bond[numBonds+extraNumBonds];
8882 memcpy(newbonds, this->bonds, numBonds*
sizeof(
Bond));
8883 memcpy(newbonds+numBonds, bonds.
begin(), extraNumBonds*
sizeof(
Bond));
8884 delete [] this->bonds;
8885 this->bonds = newbonds;
8886 numBonds += extraNumBonds;
8901 if ( extraNumAngles ) {
8902 iout <<
iINFO <<
"READ " << extraNumAngles <<
" EXTRA ANGLES\n" <<
endi;
8903 if ( anglesNormal ) {
8904 iout <<
iINFO <<
"EXTRA ANGLES ARE NORMAL HARMONIC\n" <<
endi;
8905 }
else if (
simParams->extraBondsCosAnglesSetByUser ) {
8908 iout <<
iWARN <<
"EXTRA ANGLES ARE COSINE-BASED BY DEFAULT TO MATCH PREVIOUS VERSIONS\n";
8909 iout <<
iWARN <<
"FOR NORMAL HARMONIC EXTRA ANGLES SET extraBondsCosAngles off\n" <<
endi;
8911 #ifndef MEM_OPT_VERSION 8912 Angle *newangles =
new Angle[numAngles+extraNumAngles];
8913 memcpy(newangles, this->angles, numAngles*
sizeof(
Angle));
8914 memcpy(newangles+numAngles, angles.
begin(), extraNumAngles*
sizeof(
Angle));
8915 delete [] this->angles;
8916 this->angles = newangles;
8917 numAngles += extraNumAngles;
8932 if ( extraNumDihedrals ) {
8933 iout <<
iINFO <<
"READ " << extraNumDihedrals <<
" EXTRA DIHEDRALS\n" <<
endi;
8934 #ifndef MEM_OPT_VERSION 8936 memcpy(newdihedrals, this->dihedrals, numDihedrals*
sizeof(
Dihedral));
8937 memcpy(newdihedrals+numDihedrals, dihedrals.
begin(), extraNumDihedrals*
sizeof(
Dihedral));
8938 delete [] this->dihedrals;
8939 this->dihedrals = newdihedrals;
8940 numDihedrals += extraNumDihedrals;
8955 if ( extraNumImpropers ) {
8956 iout <<
iINFO <<
"READ " << extraNumImpropers <<
" EXTRA IMPROPERS\n" <<
endi;
8957 #ifndef MEM_OPT_VERSION 8959 memcpy(newimpropers, this->impropers, numImpropers*
sizeof(
Improper));
8960 memcpy(newimpropers+numImpropers, impropers.
begin(), extraNumImpropers*
sizeof(
Improper));
8961 delete [] this->impropers;
8962 this->impropers = newimpropers;
8963 numImpropers += extraNumImpropers;
8997 PDB *initial_pdb,
char *cwd,
8998 const char *simmethod) {
9008 if (alchfile == NULL) {
9009 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, alchfile required.");
9011 strcpy(filename,
"coordinate pdb file (default)");
9014 if (alchfile->
next != NULL) {
9015 char *new_err_msg =
new char[24 + strlen(simmethod) + 26];
9016 sprintf(new_err_msg,
"Multiple definitions of %sFile in configuration file",simmethod);
9020 if ((cwd == NULL) || (alchfile->
data[0] ==
'/')) {
9021 strcpy(filename, alchfile->
data);
9024 strcpy(filename, cwd);
9025 strcat(filename, alchfile->
data);
9028 bPDB =
new PDB(filename);
9030 NAMD_die(
"Memory allocation failed in Molecule:build_fep_flags");
9034 char *new_err_msg =
new char[19 + strlen(simmethod) + 38];
9035 sprintf(new_err_msg,
"Number of atoms in %sFile PDB does not match coordinate PDB",simmethod);
9043 if (alchcol == NULL) {
9047 if (alchcol->
next != NULL) {
9048 char *new_err_msg =
new char[24 + strlen(simmethod) + 35];
9049 sprintf(new_err_msg,
"Multiple definitions of %s parameter column in config file",simmethod);
9053 if (strcasecmp(alchcol->
data,
"X") == 0) {
9056 else if (strcasecmp(alchcol->
data,
"Y") == 0) {
9059 else if (strcasecmp(alchcol->
data,
"Z") == 0) {
9062 else if (strcasecmp(alchcol->
data,
"O") == 0) {
9065 else if (strcasecmp(alchcol->
data,
"B") == 0) {
9069 NAMD_die(
"alchcol must have value of X, Y, Z, O or B");
9073 iout <<
iINFO <<
"To read " << simmethod <<
"data from file: " << filename
9075 iout <<
iINFO <<
"To read " << simmethod <<
"flag data from column: " << bcol
9079 fepAtomFlags =
new unsigned char[numAtoms];
9081 if (fepAtomFlags == NULL) {
9082 NAMD_die(
"Memory allocation failed in Molecule::build_fep_params()");
9085 double lesMassFactor = 1.0;
9087 lesMassFactor = 1.0 /
simParams->lesFactor;
9091 for (i = 0; i < numAtoms; i++) {
9095 bval = (bPDB->
atom(i))->xcoor();
9098 bval = (bPDB->
atom(i))->ycoor();
9101 bval = (bPDB->
atom(i))->zcoor();
9104 bval = (bPDB->
atom(i))->occupancy();
9107 bval = (bPDB->
atom(i))->temperaturefactor();
9113 if ( bval == (
int) bval && bval > 0 ) {
9115 NAMD_die(
"LES flag must be less than or equal to lesFactor.");
9116 fepAtomFlags[i] = (int) bval;
9117 #ifdef MEM_OPT_VERSION 9118 Real newMass = atomMassPool[eachAtomMass[i]]*lesMassFactor;
9119 eachAtomMass[i] = insert_new_mass(newMass);
9121 atoms[i].mass *= lesMassFactor;
9126 fepAtomFlags[i] = 0;
9130 fepAtomFlags[i] = 3;
9132 }
else if (bval == 1.0) {
9135 }
else if (bval == -1.0) {
9136 fepAtomFlags[i] = 2;
9138 }
else if (bval == -2.0) {
9139 fepAtomFlags[i] = 4;
9142 fepAtomFlags[i] = 0;
9144 }
else if (
simParams->pairInteractionOn) {
9145 if (bval ==
simParams->pairInteractionGroup1) {
9146 fepAtomFlags[i] = 1;
9148 }
else if (bval ==
simParams->pairInteractionGroup2) {
9149 fepAtomFlags[i] = 2;
9152 fepAtomFlags[i] = 0;
9154 }
else if (
simParams->pressureProfileAtomTypes > 1) {
9155 fepAtomFlags[i] = (int) bval;
9157 #ifdef OPENATOM_VERSION 9161 fepAtomFlags[i] = bval;
9164 fepAtomFlags[i] = 0;
9167 #endif //OPENATOM_VERSION 9171 if (alchfile != NULL) {
9195 if (ssfile == NULL) {
9196 if ( ! initial_pdb ) {
9197 NAMD_die(
"Initial PDB file unavailable, soluteScalingFile required.");
9200 strcpy(filename,
"coordinate PDB file (default)");
9203 if (ssfile->
next != NULL) {
9204 NAMD_die(
"Multiple definitions of soluteScalingFile in configuration file");
9207 if ((cwd == NULL) || (ssfile->
data[0] ==
'/')) {
9208 strcpy(filename, ssfile->
data);
9211 strcpy(filename, cwd);
9212 strcat(filename, ssfile->
data);
9215 bPDB =
new PDB(filename);
9217 NAMD_die(
"Memory allocation failed in Molecule::build_ss_flags");
9221 NAMD_die(
"Number of atoms in soluteScalingFile PDB does not match coordinate PDB");
9225 if (sscol == NULL) {
9229 if (sscol->
next != NULL) {
9230 NAMD_die(
"Multiple definitions of soluteScalingCol value in config file");
9233 if (strcasecmp(sscol->
data,
"X") == 0) {
9236 else if (strcasecmp(sscol->
data,
"Y") == 0) {
9239 else if (strcasecmp(sscol->
data,
"Z") == 0) {
9242 else if (strcasecmp(sscol->
data,
"O") == 0) {
9245 else if (strcasecmp(sscol->
data,
"B") == 0) {
9249 NAMD_die(
"soluteScalingCol must have value of X, Y, Z, O or B");
9253 iout <<
iINFO <<
"Reading solute scaling data from file: " 9254 << filename <<
"\n" <<
endi;
9255 iout <<
iINFO <<
"Reading solute scaling flags from column: " 9256 << bcol <<
"\n" <<
endi;
9258 ssAtomFlags =
new unsigned char[numAtoms];
9259 ss_index =
new int[numAtoms];
9261 if (ssAtomFlags == NULL || ss_index == NULL) {
9262 NAMD_die(
"Memory allocation failed in Molecule::build_ss_params()");
9266 for (i = 0; i < numAtoms; i++) {
9269 bval = (bPDB->
atom(i))->xcoor();
9272 bval = (bPDB->
atom(i))->ycoor();
9275 bval = (bPDB->
atom(i))->zcoor();
9278 bval = (bPDB->
atom(i))->occupancy();
9281 bval = (bPDB->
atom(i))->temperaturefactor();
9287 ss_index[num_ss] = i;
9296 if (ssfile != NULL) {
9305 int *numAtomsByLjType =
new int[LJtypecount];
9309 ss_vdw_type =
new int[LJtypecount];
9312 for (i = 0; i < LJtypecount; i++) {
9313 numAtomsByLjType[i] = 0;
9318 for (i = 0; i < num_ss; i++) {
9319 numAtomsByLjType[atoms[ss_index[i]].vdw_type]++;
9323 ss_num_vdw_params = 0;
9324 for (i = 0; i < LJtypecount; i++) {
9326 if (numAtomsByLjType[i] != 0) {
9329 ss_vdw_type[ss_num_vdw_params] = i;
9332 ss_num_vdw_params++;
9336 for (i = 0; i < num_ss; i++) {
9338 for (j = 0; j < ss_num_vdw_params; j++) {
9340 if (atoms[ss_index[i]].vdw_type == ss_vdw_type[j]) {
9343 atoms[ss_index[i]].vdw_type = LJtypecount + j;
9348 delete [] numAtomsByLjType;
9362 #ifndef MEM_OPT_VERSION 9366 FILE *alch_unpert_bond_file;
9369 if ((alch_unpert_bond_file =
Fopen(alch_fname,
"r")) == NULL) {
9370 sprintf(err_msg,
"UNABLE TO OPEN ALCH UNPERTBURBED BOND FILE %s", alch_fname);
9380 NAMD_die(
"DID NOT FIND NBOND AFTER ATOM LIST IN ALCH UNPERT PSF");
9384 sscanf(buffer,
"%d", &num_alch_unpert_Bonds);
9386 read_alch_unpert_bonds(alch_unpert_bond_file);
9395 NAMD_die(
"DID NOT FIND NTHETA AFTER BOND LIST IN ALCH UNPERT PSF");
9399 sscanf(buffer,
"%d", &num_alch_unpert_Angles);
9401 read_alch_unpert_angles(alch_unpert_bond_file);
9412 NAMD_die(
"DID NOT FIND NPHI AFTER ANGLE LIST IN ALCH UNPERT PSF");
9416 sscanf(buffer,
"%d", &num_alch_unpert_Dihedrals);
9418 read_alch_unpert_dihedrals(alch_unpert_bond_file);
9420 Fclose(alch_unpert_bond_file);
9426 suspiciousAlchBonds = 0;
9427 for (
int i = 0; i < numBonds; i++) {
9428 int part1 = fepAtomFlags[bonds[i].atom1];
9429 int part2 = fepAtomFlags[bonds[i].atom2];
9430 if ((part1 == 1 || part2 == 1 ) &&
9431 (part1 == 2 || part2 == 2 )) {
9433 suspiciousAlchBonds++;
9438 Angle *nonalchAngles;
9439 nonalchAngles =
new Angle[numAngles];
9440 int nonalchAngleCount = 0;
9441 alchDroppedAngles = 0;
9442 for (
int i = 0; i < numAngles; i++) {
9443 int part1 = fepAtomFlags[angles[i].atom1];
9444 int part2 = fepAtomFlags[angles[i].atom2];
9445 int part3 = fepAtomFlags[angles[i].atom3];
9446 if ((part1 == 1 || part2 == 1 || part3 == 1) &&
9447 (part1 == 2 || part2 == 2 || part3 == 2)) {
9449 alchDroppedAngles++;
9452 if ( angles[i].angle_type == -1 ) {
9455 "MISSING PARAMETERS FOR ANGLE %i %i %i PARTITIONS %i %i %i\n",
9456 angles[i].atom1+1, angles[i].atom2+1, angles[i].atom3+1,
9457 part1, part2, part3);
9460 nonalchAngles[nonalchAngleCount++] = angles[i];
9463 numAngles = nonalchAngleCount;
9465 angles =
new Angle[numAngles];
9466 for (
int i = 0; i < nonalchAngleCount; i++) {
9467 angles[i]=nonalchAngles[i];
9469 delete [] nonalchAngles;
9474 nonalchDihedrals =
new Dihedral[numDihedrals];
9475 int nonalchDihedralCount = 0;
9476 alchDroppedDihedrals = 0;
9477 for (
int i = 0; i < numDihedrals; i++) {
9478 int part1 = fepAtomFlags[dihedrals[i].atom1];
9479 int part2 = fepAtomFlags[dihedrals[i].atom2];
9480 int part3 = fepAtomFlags[dihedrals[i].atom3];
9481 int part4 = fepAtomFlags[dihedrals[i].atom4];
9482 if ((part1 == 1 || part2 == 1 || part3 == 1 || part4 == 1) &&
9483 (part1 == 2 || part2 == 2 || part3 == 2 || part4 == 2)) {
9485 alchDroppedDihedrals++;
9488 if ( dihedrals[i].dihedral_type == -1 ) {
9491 "MISSING PARAMETERS FOR DIHEDRAL %i %i %i %i PARTITIONS %i %i %i %i\n",
9492 dihedrals[i].atom1+1, dihedrals[i].atom2+1,
9493 dihedrals[i].atom3+1, dihedrals[i].atom4+1,
9494 part1, part2, part3, part4);
9497 nonalchDihedrals[nonalchDihedralCount++] = dihedrals[i];
9500 numDihedrals = nonalchDihedralCount;
9501 delete [] dihedrals;
9502 dihedrals =
new Dihedral[numDihedrals];
9503 for (
int i = 0; i < numDihedrals; i++) {
9504 dihedrals[i]=nonalchDihedrals[i];
9506 delete [] nonalchDihedrals;
9510 nonalchImpropers =
new Improper[numImpropers];
9511 int nonalchImproperCount = 0;
9512 alchDroppedImpropers = 0;
9513 for (
int i = 0; i < numImpropers; i++) {
9514 int part1 = fepAtomFlags[impropers[i].atom1];
9515 int part2 = fepAtomFlags[impropers[i].atom2];
9516 int part3 = fepAtomFlags[impropers[i].atom3];
9517 int part4 = fepAtomFlags[impropers[i].atom4];
9518 if ((part1 == 1 || part2 == 1 || part3 == 1 || part4 == 1) &&
9519 (part1 == 2 || part2 == 2 || part3 == 2 || part4 == 2)) {
9521 alchDroppedImpropers++;
9524 nonalchImpropers[nonalchImproperCount++] = impropers[i];
9527 numImpropers = nonalchImproperCount;
9528 delete [] impropers;
9529 impropers =
new Improper[numImpropers];
9530 for (
int i = 0; i < numImpropers; i++) {
9531 impropers[i]=nonalchImpropers[i];
9533 delete [] nonalchImpropers;
9554 if (fixedfile == NULL) {
9555 if ( ! initial_pdb )
NAMD_die(
"Initial PDB file unavailable, excludeFromPressureFile required.");
9558 if (fixedfile->
next != NULL) {
9559 NAMD_die(
"Multiple definitions of excluded pressure atoms PDB file in configuration file");
9562 if ( (cwd == NULL) || (fixedfile->
data[0] ==
'/') ) {
9563 strcpy(filename, fixedfile->
data);
9565 strcpy(filename, cwd);
9566 strcat(filename, fixedfile->
data);
9568 bPDB =
new PDB(filename);
9569 if ( bPDB == NULL ) {
9570 NAMD_die(
"Memory allocation failed in Molecule::build_exPressure_atoms");
9574 NAMD_die(
"Number of atoms in excludedPressure atoms PDB doesn't match coordinate PDB");
9583 if (fixedcol == NULL) {
9586 if (fixedcol->
next != NULL) {
9587 NAMD_die(
"Multiple definitions of excludedPressure atoms column in config file");
9590 if (strcasecmp(fixedcol->
data,
"X") == 0) {
9592 }
else if (strcasecmp(fixedcol->
data,
"Y") == 0) {
9594 }
else if (strcasecmp(fixedcol->
data,
"Z") == 0) {
9596 }
else if (strcasecmp(fixedcol->
data,
"O") == 0) {
9598 }
else if (strcasecmp(fixedcol->
data,
"B") == 0) {
9601 NAMD_die(
"excludedPressureFileCol must have value of X, Y, Z, O, or B");
9606 exPressureAtomFlags =
new int32[numAtoms];
9608 if (exPressureAtomFlags == NULL) {
9609 NAMD_die(
"memory allocation failed in Molecule::build_fixed_atoms()");
9612 numExPressureAtoms = 0;
9615 for (i=0; i<numAtoms; i++) {
9618 case 1: bval = (bPDB->
atom(i))->xcoor();
break;
9619 case 2: bval = (bPDB->
atom(i))->ycoor();
break;
9620 case 3: bval = (bPDB->
atom(i))->zcoor();
break;
9621 case 4: bval = (bPDB->
atom(i))->occupancy();
break;
9622 case 5: bval = (bPDB->
atom(i))->temperaturefactor();
break;
9627 exPressureAtomFlags[i] = 1;
9628 numExPressureAtoms++;
9630 exPressureAtomFlags[i] = 0;
9633 if (fixedfile != NULL)
9636 iout <<
iINFO <<
"Got " << numExPressureAtoms <<
" excluded pressure atoms." 9646 return ((atoms[anum].status &
DrudeAtom) != 0);
9656 return ((atoms[anum].status &
OxygenAtom) != 0);
9661 return (hydrogenGroup[atoms[anum].hydrogenList].isGP);
9666 return (hydrogenGroup[atoms[anum].hydrogenList].waterVal == 2);
9671 return (hydrogenGroup[atoms[anum].hydrogenList].atomsInGroup);
9678 return atoms[anum].partner;
9682 if ( n != numAtoms )
9683 NAMD_die(
"Incorrect number of atoms in Molecule::reloadCharges().");
9685 #ifdef MEM_OPT_VERSION 9686 delete [] atomChargePool;
9687 vector<Real> tmpCharges;
9688 for(
int i=0; i<numAtoms; i++){
9692 for(
int j=0; j<tmpCharges.size();j++){
9693 if(tmpCharges[j] == charge[i]){
9699 tmpCharges.push_back(charge[i]);
9700 foundIdx = tmpCharges.size()-1;
9702 eachAtomCharge[i] = (
Index)foundIdx;
9704 chargePoolSize = tmpCharges.size();
9705 atomChargePool =
new Real[chargePoolSize];
9706 for(
int i=0; i<chargePoolSize; i++)
9707 atomChargePool[i] = tmpCharges[i];
9709 for(
int i=0; i<n; ++i ) atoms[i].charge = charge[i];
9720 PDB dcdSelectionInputFilePdb(dcdSelectionInputFile);
9721 int numDcdSelectionAtoms = dcdSelectionInputFilePdb.num_atoms();
9723 if (numDcdSelectionAtoms < 1)
9724 NAMD_die(
"No atoms found in dcdSelectionInputFile\n");
9729 if (numDcdSelectionAtoms != numAtoms)
9730 NAMD_die(
"The number of atoms in dcdSelectionInputFile must be equal to the total number of atoms in the structure!");
9731 int bitmask = 1 << index;
9732 for (
int i=0; i<numAtoms; i++) {
9733 #ifdef MEM_OPT_VERSION 9734 PDBCoreData *atom = dcdSelectionInputFilePdb.atom(i);
9736 PDBAtom *atom = dcdSelectionInputFilePdb.atom(i);
9741 atoms[i].flags.dcdSelection |= bitmask;
9744 dcdSelectionParams[index].tag = index+1;
9750 dcdSelectionParams[index].dcdSelectionIndexReverse.resize(numAtoms);
9752 int bitmask = 1 << index;
9754 for (
int i=0; i<numAtoms; i++) {
9755 if(atoms[i].flags.dcdSelection & bitmask)
9757 dcdSelectionParams[index].dcdSelectionIndex.push_back(i);
9758 dcdSelectionParams[index].dcdSelectionIndexReverse[i] = count++;
9762 dcdSelectionParams[index].dcdSelectionIndexReverse[i] = -1;
9765 dcdSelectionParams[index].size = count;
9770 dcdSelectionParams[index].frequency = freq;
9776 if(
auto keyIt=dcdSelectionKeyMap.find(std::string(keystr)); keyIt!=dcdSelectionKeyMap.end())
9778 return(keyIt->second);
9786 uint16 key=find_dcd_selection_index(keystr);
9793 key = dcdSelectionKeyMap.size();
9797 sprintf(errmsg,
"Key %s beyond Max (16) DCD Selections\n", keystr);
9800 dcdSelectionKeyMap[std::string(keystr)] = key;
9816 char *keystr =
new char[256];
9817 char *valstr =
new char[256];
9820 int dcdSelectionCounter[3]={0,0,0};
9821 StringList *dcdSelectionInputFileList = config->
find(
"dcdselectioninputfile");
9822 for(;dcdSelectionInputFileList; dcdSelectionInputFileList = dcdSelectionInputFileList->
next)
9824 sscanf(dcdSelectionInputFileList->
data,
"%80s %255s",keystr,valstr);
9825 key = find_or_create_dcd_selection_index(keystr);
9826 dcdSelectionCounter[0]++;
9827 build_dcd_selection_list_pdb(key, valstr);
9828 iout <<
iINFO <<
"Dcdselection tag "<< keystr <<
" Input file " << valstr <<
" index "<< key <<
"\n";
9830 StringList *dcdSelectionFileList = config->
find(
"dcdselectionfile");
9831 for(;dcdSelectionFileList; dcdSelectionFileList = dcdSelectionFileList->
next)
9834 sscanf(dcdSelectionFileList->
data,
"%80s %255s",keystr,valstr);
9835 key = find_dcd_selection_index(keystr);
9838 add_dcd_selection_file(key, valstr);
9839 dcdSelectionCounter[1]++;
9843 sprintf(errmsg,
"Keystring %s has no dcd input file match\n", keystr);
9846 iout <<
iINFO <<
"Dcdselection tag "<< keystr <<
" Output file " << valstr<<
" index "<< key <<
"\n";
9848 StringList *dcdSelectionFreqList = config->
find(
"dcdselectionfreq");
9849 for(;dcdSelectionFreqList; dcdSelectionFreqList = dcdSelectionFreqList->
next)
9852 sscanf(dcdSelectionFreqList->
data,
"%s %d",keystr,&freq);
9853 key = find_dcd_selection_index(keystr);
9856 add_dcd_selection_freq(key, freq);
9857 dcdSelectionCounter[2]++;
9861 sprintf(errmsg,
"Keystring %s has no dcd input file match\n", keystr);
9864 iout <<
iINFO <<
"Dcd Selection tag "<< keystr <<
" freq " << freq<<
" index "<< key <<
"\n";
9869 if((dcdSelectionCounter[0]!=dcdSelectionCounter[1] && CmiNumPartitions()==1) || dcdSelectionCounter[0]!=dcdSelectionCounter[2] || dcdSelectionCounter[0]>15)
9872 sprintf(errmsg,
"Invalid DCD Selection configuration\n");
9878 #ifndef MEM_OPT_VERSION 9885 void Molecule::build_atom_status(
void) {
9888 int numDrudeWaters = 0;
9891 numLonepairs = numZeroMassAtoms;
9893 iout <<
iWARN <<
"CORRECTION OF ZERO MASS ATOMS TURNED OFF " 9894 "BECAUSE LONE PAIRS ARE USED\n" <<
endi;
9898 if (is_lonepairs_psf && numLonepairs != numLphosts) {
9899 NAMD_die(
"must have same number of LP hosts as lone pairs");
9901 }
else if (numZeroMassAtoms) {
9902 for (i=0; i < numAtoms; i++) {
9903 if ( atoms[i].mass < 0.001 ) atoms[i].mass = 0.001;
9906 iout <<
iWARN <<
"FOUND " << numZeroMassAtoms <<
9907 " ATOMS WITH ZERO OR NEGATIVE MASSES! CHANGED TO 0.001\n" <<
endi;
9912 hydrogenGroup.resize(numAtoms);
9914 for (i=0; i < numAtoms; i++) {
9915 atoms[i].partner = (-1);
9925 int hhbondcount = 0;
9926 for (i=0; i < numRealBonds; i++) {
9927 a1 = bonds[i].atom1;
9928 a2 = bonds[i].atom2;
9929 if (is_hydrogen(a1) && is_hydrogen(a2)) {
9932 atoms[a1].partner = a2;
9933 atoms[a2].partner = a1;
9941 if ( hhbondcount && ! CkMyPe() ) {
9942 iout <<
iWARN <<
"Found " << hhbondcount <<
" H-H bonds.\n" <<
endi;
9947 for (i=0; i < numRealBonds; i++) {
9948 a1 = bonds[i].atom1;
9949 a2 = bonds[i].atom2;
9950 if (is_hydrogen(a1)) {
9951 if (is_hydrogen(a2))
continue;
9952 atoms[a1].partner = a2;
9958 if (is_oxygen(a2)) hg[a2].waterVal++;
9960 if (is_hydrogen(a2)) {
9961 atoms[a2].partner = a1;
9967 if (is_oxygen(a1)) hg[a1].waterVal++;
9973 atoms[a1].partner = a2;
9980 atoms[a2].partner = a1;
9988 else if ( is_lonepairs_psf || is_drude_psf ) {
9989 if (is_lp(a1) || is_drude(a1)) {
9990 if (is_hydrogen(a2) || is_lp(a2) || is_drude(a2)) {
9992 sprintf(msg,
"%s particle %d is bonded to non-parent atom %d",
9993 (is_lp(a1) ?
"Lone pair" :
"Drude"), a1+1, a2+1);
9996 atoms[a1].partner = a2;
10002 else if (is_lp(a2) || is_drude(a2)) {
10003 if (is_hydrogen(a1) || is_lp(a1) || is_drude(a1)) {
10005 sprintf(msg,
"%s particle %d is bonded to non-parent atom %d",
10006 (is_lp(a2) ?
"Lone pair" :
"Drude"), a2+1, a1+1);
10009 atoms[a2].partner = a1;
10021 for(i=0; i<numAtoms; i++) {
10022 if ( ! hg[hg[i].GPID].isGP ) {
10024 sprintf(msg,
"child atom %d bonded only to child H atoms",i+1);
10027 if ( hg[i].isGP && is_hydrogen(i) ) {
10028 if ( hg[i].GPID == i )
continue;
10030 if ( is_hydrogen(hg[i].GPID) && hg[hg[i].GPID].GPID != i ) {
10032 sprintf(msg,
"H atom %d bonded only to child H atoms",i+1);
10038 if ( hg[i].atomsInGroup != 2 ) {
10040 sprintf(msg,
"H atom %d bonded to multiple H atoms",i+1);
10045 if ( hGPcount && ! CkMyPe() ) {
10046 iout <<
iWARN <<
"Found " << hGPcount <<
" H-H molecules.\n" <<
endi;
10050 for (i=0; i<numAtoms; ++i) {
10051 if ( hg[i].isGP ) hg[i].
GPID = i;
10059 for (i=0; i<numLphosts; ++i) {
10060 int a1 = lphosts[i].atom1;
10061 int a2 = lphosts[i].atom2;
10062 int a3 = lphosts[i].atom3;
10063 int a4 = lphosts[i].atom4;
10064 int m1 = hg[a1].
MPID;
10065 while ( hg[m1].MPID != m1 ) m1 = hg[m1].
MPID;
10066 int m2 = hg[a2].
MPID;
10067 while ( hg[m2].MPID != m2 ) m2 = hg[m2].
MPID;
10068 int m3 = hg[a3].
MPID;
10069 while ( hg[m3].MPID != m3 ) m3 = hg[m3].
MPID;
10070 int m4 = hg[a4].
MPID;
10071 while ( hg[m4].MPID != m4 ) m4 = hg[m4].
MPID;
10073 if ( m2 < mp ) mp = m2;
10074 if ( m3 < mp ) mp = m3;
10075 if ( m4 < mp ) mp = m4;
10083 for (i=0; i<numAtoms; ++i) {
10084 int mp = hg[i].
MPID;
10085 if ( hg[mp].MPID != mp ) {
10090 if ( allok )
break;
10092 for (i=0; i<numAtoms; ++i) {
10096 for (i=0; i<numAtoms; ++i) {
10102 for (i=0; i<numAtoms; i++) {
10113 numHydrogenGroups = 0;
10114 maxHydrogenGroupSize = 0;
10115 numMigrationGroups = 0;
10116 maxMigrationGroupSize = 0;
10117 for(i=0; i<numAtoms; i++)
10120 ++numMigrationGroups;
10122 if ( mgs > maxMigrationGroupSize ) maxMigrationGroupSize = mgs;
10125 ++numHydrogenGroups;
10127 if ( hgs > maxHydrogenGroupSize ) maxHydrogenGroupSize = hgs;
10131 hydrogenGroup.sort();
10136 for(i=0; i<numAtoms; ++i, --hgs) {
10138 if ( hg[i].isGP ) {
10140 parentid = hg[i].
atomID;
10143 sprintf(buff,
"Atom %d has bad hydrogen group size. " 10144 "Check for duplicate bonds.", parentid+1);
10148 if ( hg[i].isGP ) {
10150 sprintf(buff,
"Atom %d has bad hydrogen group size. " 10151 "Check for duplicate bonds.", parentid+1);
10159 for(i=0; i<numAtoms; ++i, --mgs) {
10161 if ( hg[i].isMP ) {
10163 parentid = hg[i].
atomID;
10166 sprintf(buff,
"Atom %d has bad migration group size.", parentid+1);
10170 if ( hg[i].isMP ) {
10172 sprintf(buff,
"Atom %d has bad migration group size.", parentid+1);
10180 for(i=0; i<numAtoms; i++) {
10181 atoms[hydrogenGroup[i].atomID].hydrogenList = i;
10187 for (i = 0; i < numAtoms; i++) {
10188 if (is_water(hg[i].atomID) && hg[i].isGP) {
10190 || ! is_drude(hg[i+1].atomID)
10191 || ! is_lp(hg[i+2].atomID)
10192 || ! is_hydrogen(hg[i+3].atomID)
10193 || ! is_hydrogen(hg[i+4].atomID) ) {
10195 sprintf(msg,
"Drude water molecule from HydrogenGroup i=%d " 10196 "starting at atom %d is not sorted\n", i, hg[i].atomID+1);
10203 else if (is_drude(hg[i].atomID)) {
10204 if (i < 1 || hg[i-1].atomID != hg[i].GPID) {
10206 sprintf(msg,
"Drude particle from HydrogenGroup i=%d must " 10207 "immediately follow its parent atom %d\n", i, hg[i].GPID+1);
10212 else if (is_lp(hg[i].atomID)) {
10214 sprintf(msg,
"Drude lonepair from HydrogenGroup i=%d " 10215 "at particle %d is NOT from water - unsupported\n",
10216 i, hg[i].atomID+1);
10226 for(i=0; i<numAtoms; i++)
10227 iout << i <<
" atomID=" << hydrogenGroup[i].atomID
10228 <<
" isGP=" << hydrogenGroup[i].isGP
10229 <<
" parent=" << hydrogenGroup[i].GPID
10230 <<
" #" << hydrogenGroup[i].atomsInGroup
10231 <<
" waterVal=" << hydrogenGroup[i].waterVal
10232 <<
" partner=" << atoms[i].partner
10233 <<
" hydrogenList=" << atoms[i].hydrogenList
10244 delete [] rigidBondLengths;
10245 rigidBondLengths =
new Real[numAtoms];
10246 if ( ! rigidBondLengths ) {
10247 NAMD_die(
"Memory allocation failed in Molecule::build_atom_status()\n");
10249 for (i=0; i<numAtoms; ++i) rigidBondLengths[i] = 0;
10254 for (i=0; i < numRealBonds; i++) {
10255 a1 = bonds[i].atom1;
10256 a2 = bonds[i].atom2;
10259 if (is_hydrogen(a2)) {
int tmp = a1; a1 = a2; a2 = tmp; }
10260 if (is_hydrogen(a1)) {
10261 if ( is_hydrogen(a2) ) {
10262 if ( ! is_water(a2) ) {
10263 rigidBondLengths[a1] = ( mode ==
RIGID_ALL ? x0 : 0. );
10264 rigidBondLengths[a2] = ( mode ==
RIGID_ALL ? x0 : 0. );
10266 }
else if ( is_water(a2) || mode ==
RIGID_ALL ) {
10267 rigidBondLengths[a1] = x0;
10268 if (is_water(a2)) r_oh = rigidBondLengths[a1];
10270 rigidBondLengths[a1] = 0.;
10275 if (is_lp(a2)) {
int tmp = a1; a1 = a2; a2 = tmp; }
10277 if (! is_water(a2) ) {
10281 sprintf(err_msg,
"ILLEGAL LONE PAIR AT INDEX %i\n" 10282 "LONE PAIRS ARE CURRENTLY ALLOWED ONLY ON WATER MOLECULES\n",
10286 rigidBondLengths[a1] = x0;
10295 int tmp = a1; a1 = a2; a2 = tmp;
10298 if (is_water(a2)) {
10305 sprintf(msg,
"ILLEGAL LONE PAIR AT INDEX %d\n" 10306 "LONE PAIRS ARE CURRENTLY ALLOWED ONLY ON WATER MOLECULES\n",
10316 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
10317 for( ; h_i != h_e; ++h_i ) {
10318 if ( h_i->
isGP ) rigidBondLengths[h_i->
atomID] = 0.;
10322 for (i=0; i < numAngles; i++) {
10323 a2 = angles[i].atom2;
10324 if ( ! is_water(a2) )
continue;
10325 if ( ! is_oxygen(a2) )
continue;
10326 a1 = angles[i].atom1;
10327 if ( ! is_hydrogen(a1) )
continue;
10328 a3 = angles[i].atom3;
10329 if ( ! is_hydrogen(a3) )
continue;
10330 if (is_lp(a2) || is_lp(a1) || is_lp(a3) ||
10331 is_drude(a2) || is_drude(a1) || is_drude(a3))
continue;
10332 if ( rigidBondLengths[a1] != rigidBondLengths[a3] ) {
10333 if (rigidBondLengths[a1] >0.3 && rigidBondLengths[a3] >0.3) {
10334 printf(
"length1: %f length2: %f\n", rigidBondLengths[a1], rigidBondLengths[a3]);
10336 NAMD_die(
"Asymmetric water molecule found??? This can't be right.\n");
10341 rigidBondLengths[a2] = 2. * rigidBondLengths[a1] * sin(0.5*t0);
10342 r_hh = rigidBondLengths[a2];
10346 int numBondWaters = 0;
10347 int numFailedWaters = 0;
10349 for (i=0; i < numRealBonds; i++) {
10350 a1 = bonds[i].atom1;
10351 a2 = bonds[i].atom2;
10352 if ( ! is_hydrogen(a1) )
continue;
10353 if ( ! is_hydrogen(a2) )
continue;
10354 int ma1 = get_mother_atom(a1);
10355 int ma2 = get_mother_atom(a2);
10356 if ( ma1 != ma2 )
continue;
10357 if ( ! is_water(ma1) )
continue;
10358 if ( rigidBondLengths[ma1] != 0. )
continue;
10361 rigidBondLengths[ma1] = x0;
10368 if (r_oh < 0.0 || r_hh < 0.0) {
10370 NAMD_die(
"Failed to find water bond lengths\n");
10372 r_ohc = sqrt(r_oh * r_oh - 0.25 * r_hh * r_hh);
10376 h_i = hydrogenGroup.begin(); h_e = hydrogenGroup.end();
10377 for( ; h_i != h_e; ++h_i ) {
10379 rigidBondLengths[h_i->
atomID] == 0. ) {
10380 if ( h_i + 1 == h_e || h_i + 2 == h_e ||
10381 h_i[1].isGP || h_i[2].isGP || h_i->
atomsInGroup != 3 ) {
10382 NAMD_die(
"Abnormal water detected.");
10384 if ( CkNumNodes() > 1 ) {
10385 NAMD_die(
"Unable to determine H-H distance for rigid water because structure has neither H-O-H angle nor H-H bond.");
10391 get_atomtype(btmp.
atom1),
10392 get_atomtype(btmp.
atom2),
10398 rigidBondLengths[h_i->
atomID] = x0;
10405 if ( numBondWaters + numFailedWaters ) {
10406 iout <<
iWARN <<
"Missing angles for " <<
10407 ( numBondWaters + numFailedWaters ) <<
" waters.\n" <<
endi;
10409 if ( numBondWaters ) {
10410 iout <<
iWARN <<
"Obtained H-H distance from bond parameters for " <<
10411 numBondWaters <<
" waters.\n" <<
endi;
10412 iout <<
iWARN <<
"This would not be possible in a multi-process run.\n" <<
endi;
10414 if ( numFailedWaters ) {
10415 iout <<
iERROR <<
"Failed to obtain H-H distance from angles or bonds for " <<
10416 numFailedWaters <<
" waters.\n" <<
endi;
10424 for (i=0; i<numAtoms; ++i) rigidBondLengths[i] *= -1;
10426 for (i=0; i<numAtoms; ++i) {
10427 if ( ! is_water(i) ) rigidBondLengths[i] *= -1;
10433 for (i=0; i<numAtoms; ++i) {
10434 if ( rigidBondLengths[i] > 0. ) ++numRigidBonds;
10455 BigReal LJAvgA, LJAvgB, LJAvgA1, LJAvgB1, LJAvgA2, LJAvgB2;
10456 int numLJsites, numLJsites1, numLJsites2;
10468 const int ss_dim = soluteScalingOn ? ss_num_vdw_params : 0;
10476 if (soluteScalingOn && (ss_vdw_type == NULL)) {
10477 NAMD_bug(
"Solute scaling is used but ss_vdw_type is NULL. " 10478 "It is likely that compute_LJcorrection() is called " 10479 "before build_ss_flags().");
10481 Real A, B, A14, B14;
10482 Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
10483 Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
10484 Real *ATable =
new Real[LJtypecount*LJtypecount];
10485 Real *BTable =
new Real[LJtypecount*LJtypecount];
10486 int useGeom =
simParams->vdwGeometricSigma;
10488 for (
int i = 0; i < LJtypecount; i++) {
10489 for (
int j = 0; j < LJtypecount; j++) {
10490 const int i_type = soluteScalingOn ? ((i >= table_dim_org) ? ss_vdw_type[i-table_dim_org] : i) : i;
10491 const int j_type = soluteScalingOn ? ((j >= table_dim_org) ? ss_vdw_type[j-table_dim_org] : j) : j;
10493 ATable[i*LJtypecount + j] = A;
10494 BTable[i*LJtypecount + j] = B;
10497 params->
get_vdw_params(&sigma_i,&epsilon_i,&sigma_i14,&epsilon_i14,i_type);
10498 params->
get_vdw_params(&sigma_j,&epsilon_j,&sigma_j14,&epsilon_j14,j_type);
10500 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10501 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10502 sigma_ij *= sigma_ij*sigma_ij;
10503 sigma_ij *= sigma_ij;
10504 ATable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10505 BTable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij;
10510 int *numAtomsByLjType =
new int[LJtypecount];
10511 for (
int i=0; i < LJtypecount; i++) {numAtomsByLjType[i]=0;}
10512 for (
int i=0; i < numAtoms; i++) {numAtomsByLjType[atoms[i].vdw_type]++;}
10519 for (
int i=0; i < LJtypecount; i++) {
10520 for (
int j=0; j < LJtypecount; j++) {
10521 A = ATable[i*LJtypecount + j];
10522 B = BTable[i*LJtypecount + j];
10523 if (!A && !B)
continue;
10524 npairs = (numAtomsByLjType[i] - int(i==j))*
BigReal(numAtomsByLjType[j]);
10525 sumOfAs += npairs*A;
10526 sumOfBs += npairs*B;
10528 if (i==j) numLJsites += numAtomsByLjType[i];
10531 delete [] numAtomsByLjType;
10535 LJAvgA = sumOfAs / count;
10536 LJAvgB = sumOfBs / count;
10557 numLJsites1 = numLJsites2 = numLJsites;
10558 int alch_counter = 0;
10559 for (
int i=0; i < numAtoms; ++i) {
10561 if (get_fep_type(i) == 2 || get_fep_type(i) == 4) alchFlagi = -1;
10562 if (get_fep_type(i) == 1 || get_fep_type(i) == 3) alchFlagi = 1;
10564 &A, &B, &A14, &B14)) {
10568 &epsilon_i14, atoms[i].vdw_type);
10570 useGeom ? sqrt(sigma_i*sigma_i) : 0.5*(sigma_i+sigma_i);
10571 BigReal epsilon_ii = sqrt(epsilon_i*epsilon_i);
10573 sigma_ii *= sigma_ii*sigma_ii;
10574 sigma_ii *= sigma_ii;
10575 A = 4.0*sigma_ii*epsilon_ii*sigma_ii;
10576 B = 4.0*sigma_ii*epsilon_ii;
10579 if (alchFlagi == 1) numLJsites2--;
10580 else if (alchFlagi == -1) numLJsites1--;
10582 for (
int j=i+1; j < numAtoms; ++j) {
10584 if (get_fep_type(j) == 2 || get_fep_type(j) == 4) alchFlagj = -1;
10585 if (get_fep_type(j) == 1 || get_fep_type(j) == 3) alchFlagj = 1;
10586 int alchFlagSum = alchFlagi + alchFlagj;
10589 if (alchFlagi == 0 && alchFlagj == 0)
continue;
10592 &A, &B, &A14, &B14)) {
10596 &epsilon_i14, atoms[i].vdw_type);
10598 &epsilon_j14, atoms[j].vdw_type);
10600 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10601 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10603 sigma_ij *= sigma_ij*sigma_ij;
10604 sigma_ij *= sigma_ij;
10605 A = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10606 B = 4.0*sigma_ij*epsilon_ij;
10608 if (!A && !B)
continue;
10613 if ( alchFlagSum > 0 ){
10618 else if ( alchFlagSum < 0 ){
10634 if ( alchFlagi == 1 || alchFlagi == -1 ) alch_counter++;
10635 if ( alch_counter == (numFepInitial + numFepFinal) )
break;
10637 LJAvgA = sumOfAs / count;
10638 LJAvgB = sumOfBs / count;
10640 LJAvgA1 = sumOfAs1 / count1;
10641 LJAvgB1 = sumOfBs1 / count1;
10647 LJAvgA2 = sumOfAs2 / count2;
10648 LJAvgB2 = sumOfBs2 / count2;
10653 if ( ! CkMyPe() ) {
10654 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 10655 <<
"ENERGY AND PRESSURE\n" <<
endi;
10656 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A0 AND B0 COEFFICIENTS " 10657 << LJAvgA2 <<
" AND " << LJAvgB2 <<
"\n" <<
endi;
10658 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A1 AND B1 COEFFICIENTS " 10659 << LJAvgA1 <<
" AND " << LJAvgB1 <<
"\n" <<
endi;
10661 numLJsites = (numLJsites1 + numLJsites2 - numLJsites);
10662 LJAvgA1 *=
BigReal(numLJsites1)*numLJsites1;
10663 LJAvgB1 *=
BigReal(numLJsites1)*numLJsites1;
10664 LJAvgA2 *=
BigReal(numLJsites2)*numLJsites2;
10665 LJAvgB2 *=
BigReal(numLJsites2)*numLJsites2;
10668 LJAvgA1 = LJAvgB1 = LJAvgA2 = LJAvgB2 = 0;
10670 if ( ! CkMyPe() ) {
10671 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 10672 <<
"ENERGY AND PRESSURE\n" <<
endi;
10673 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A AND B COEFFICIENTS " 10674 << LJAvgA <<
" AND " << LJAvgB <<
"\n" <<
endi;
10677 LJAvgA *=
BigReal(numLJsites)*numLJsites;
10678 LJAvgB *=
BigReal(numLJsites)*numLJsites;
10687 BigReal rswitch2 = rswitch*rswitch;
10688 BigReal rswitch3 = rswitch*rswitch2;
10689 BigReal rswitch4 = rswitch2*rswitch2;
10690 BigReal rswitch5 = rswitch2*rswitch3;
10691 BigReal rswitch6 = rswitch3*rswitch3;
10717 BigReal int_U_gofr_A, int_rF_gofr_A, int_U_gofr_B, int_rF_gofr_B;
10736 BigReal rsum3 = (rcut + rswitch)*(rcut + rswitch)*(rcut + rswitch);
10737 int_U_gofr_A = int_rF_gofr_A = (16*
PI*(3*rcut4 + 9*rcut3*rswitch
10738 + 11*rcut2*rswitch2 + 9*rcut*rswitch3
10740 / (315*rcut5*rswitch5*rsum3));
10741 int_U_gofr_B = int_rF_gofr_B = -16*
PI / (3*rsum3);
10768 BigReal lnr = log(rcut/rswitch);
10792 int_U_gofr_A = int_rF_gofr_A =\
10793 16*
PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
10794 int_U_gofr_B = int_rF_gofr_B = -4*
PI*lnr / (rcut3 - rswitch3);
10831 int_rF_gofr_A = 8*
PI / (9*rcut9);
10832 int_rF_gofr_B = -4*
PI / (3*rcut3);
10833 int_U_gofr_A = 2*
PI / (9*rcut9);
10834 int_U_gofr_B = -2*
PI / (3*rcut3);
10837 tail_corr_virial = int_rF_gofr_A*LJAvgA + int_rF_gofr_B*LJAvgB;
10838 tail_corr_ener = int_U_gofr_A*LJAvgA + int_U_gofr_B*LJAvgB;
10840 tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1 -
10842 tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1 -
10844 tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2 -
10846 tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2 -
10865 BigReal LJAvgA, LJAvgB, LJAvgA1, LJAvgB1, LJAvgA2, LJAvgB2;
10866 int numLJsites, numLJsites1, numLJsites2;
10878 const int ss_dim = soluteScalingOn ? ss_num_vdw_params : 0;
10886 if (soluteScalingOn && (ss_vdw_type == NULL)) {
10887 NAMD_bug(
"Solute scaling is used but ss_vdw_type is NULL. " 10888 "It is likely that compute_LJcorrection_alternative() is called " 10889 "before build_ss_flags().");
10891 Real A, B, A14, B14;
10892 Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
10893 Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
10894 Real *ATable =
new Real[LJtypecount*LJtypecount];
10895 Real *BTable =
new Real[LJtypecount*LJtypecount];
10896 int useGeom =
simParams->vdwGeometricSigma;
10898 for (
int i = 0; i < LJtypecount; i++) {
10899 for (
int j = 0; j < LJtypecount; j++) {
10900 const int i_type = soluteScalingOn ? ((i >= table_dim_org) ? ss_vdw_type[i-table_dim_org] : i) : i;
10901 const int j_type = soluteScalingOn ? ((j >= table_dim_org) ? ss_vdw_type[j-table_dim_org] : j) : j;
10903 ATable[i*LJtypecount + j] = A;
10904 BTable[i*LJtypecount + j] = B;
10907 params->
get_vdw_params(&sigma_i,&epsilon_i,&sigma_i14,&epsilon_i14,i_type);
10908 params->
get_vdw_params(&sigma_j,&epsilon_j,&sigma_j14,&epsilon_j14,j_type);
10910 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10911 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10912 sigma_ij *= sigma_ij*sigma_ij;
10913 sigma_ij *= sigma_ij;
10914 ATable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10915 BTable[i*LJtypecount + j] = 4.0*sigma_ij*epsilon_ij;
10920 int *numAtomsByLjType =
new int[LJtypecount];
10921 for (
int i=0; i < LJtypecount; i++) {numAtomsByLjType[i]=0;}
10922 for (
int i=0; i < numAtoms; i++) {numAtomsByLjType[atoms[i].vdw_type]++;}
10927 for (
int i=0; i < LJtypecount; i++) {
10928 for (
int j=0; j < LJtypecount; j++) {
10929 A = ATable[i*LJtypecount + j];
10930 B = BTable[i*LJtypecount + j];
10931 npairs =
BigReal(numAtomsByLjType[i])*
BigReal(numAtomsByLjType[j]);
10932 sumOfAs += npairs*A;
10933 sumOfBs += npairs*B;
10936 delete [] numAtomsByLjType;
10960 int alch_counter = 0;
10961 for (
int i=0; i < numAtoms; ++i) {
10963 if (get_fep_type(i) == 2 || get_fep_type(i) == 4) alchFlagi = -1;
10964 if (get_fep_type(i) == 1 || get_fep_type(i) == 3) alchFlagi = 1;
10966 for (
int j=i; j < numAtoms; ++j) {
10968 if (get_fep_type(j) == 2 || get_fep_type(j) == 4) alchFlagj = -1;
10969 if (get_fep_type(j) == 1 || get_fep_type(j) == 3) alchFlagj = 1;
10970 int alchFlagSum = alchFlagi + alchFlagj;
10973 if (alchFlagi == 0 && alchFlagj == 0)
continue;
10976 &A, &B, &A14, &B14)) {
10980 &epsilon_i14, atoms[i].vdw_type);
10982 &epsilon_j14, atoms[j].vdw_type);
10984 useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
10985 BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
10987 sigma_ij *= sigma_ij*sigma_ij;
10988 sigma_ij *= sigma_ij;
10989 A = 4.0*sigma_ij*epsilon_ij*sigma_ij;
10990 B = 4.0*sigma_ij*epsilon_ij;
10995 if ( alchFlagSum > 0 ){
10998 }
else if ( alchFlagSum < 0 ){
11010 if ( alchFlagi == 1 || alchFlagi == -1 ) alch_counter++;
11011 if ( alch_counter == (numFepInitial + numFepFinal) )
break;
11016 LJAvgA1 = sumOfAs1;
11017 LJAvgB1 = sumOfBs1;
11018 LJAvgA2 = sumOfAs2;
11019 LJAvgB2 = sumOfBs2;
11021 if ( ! CkMyPe() ) {
11022 BigReal inv_sizeSq = 1.0 / (numAtoms * numAtoms);
11023 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 11024 <<
"ENERGY AND PRESSURE\n" <<
endi;
11025 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A0 AND B0 COEFFICIENTS " 11026 << LJAvgA2*inv_sizeSq <<
" AND " << LJAvgB2*inv_sizeSq <<
"\n" <<
endi;
11027 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A1 AND B1 COEFFICIENTS " 11028 << LJAvgA1*inv_sizeSq <<
" AND " << LJAvgB1*inv_sizeSq <<
"\n" <<
endi;
11031 LJAvgA1 = LJAvgB1 = LJAvgA2 = LJAvgB2 = 0;
11033 if ( ! CkMyPe() ) {
11034 BigReal inv_sizeSq = 1.0 / (numAtoms * numAtoms);
11035 iout <<
iINFO <<
"LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO " 11036 <<
"ENERGY AND PRESSURE\n" <<
endi;
11037 iout <<
iINFO <<
"LONG-RANGE LJ: AVERAGE A AND B COEFFICIENTS " 11038 << LJAvgA*inv_sizeSq <<
" AND " << LJAvgB*inv_sizeSq <<
"\n" <<
endi;
11049 BigReal rswitch2 = rswitch*rswitch;
11050 BigReal rswitch3 = rswitch*rswitch2;
11051 BigReal rswitch4 = rswitch2*rswitch2;
11052 BigReal rswitch5 = rswitch2*rswitch3;
11053 BigReal rswitch6 = rswitch3*rswitch3;
11079 BigReal int_U_gofr_A, int_rF_gofr_A, int_U_gofr_B, int_rF_gofr_B;
11098 BigReal rsum3 = (rcut + rswitch)*(rcut + rswitch)*(rcut + rswitch);
11099 int_U_gofr_A = int_rF_gofr_A = (16*
PI*(3*rcut4 + 9*rcut3*rswitch
11100 + 11*rcut2*rswitch2 + 9*rcut*rswitch3
11102 / (315*rcut5*rswitch5*rsum3));
11103 int_U_gofr_B = int_rF_gofr_B = -16*
PI / (3*rsum3);
11130 BigReal lnr = log(rcut/rswitch);
11154 int_U_gofr_A = int_rF_gofr_A =\
11155 16*
PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
11156 int_U_gofr_B = int_rF_gofr_B = -4*
PI*lnr / (rcut3 - rswitch3);
11193 int_rF_gofr_A = 8*
PI / (9*rcut9);
11194 int_rF_gofr_B = -4*
PI / (3*rcut3);
11195 int_U_gofr_A = 2*
PI / (9*rcut9);
11196 int_U_gofr_B = -2*
PI / (3*rcut3);
11199 tail_corr_virial = int_rF_gofr_A*LJAvgA + int_rF_gofr_B*LJAvgB;
11200 tail_corr_ener = int_U_gofr_A*LJAvgA + int_U_gofr_B*LJAvgB;
11202 tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1 -
11204 tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1 -
11206 tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2 -
11208 tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2 -
11220 const BigReal corr = (doti ? 0.0 : tail_corr_ener);
11221 return scl*(corr + vdw_lambda_1*tail_corr_dUdl_1 +
11222 vdw_lambda_2*tail_corr_dUdl_2);
11225 return scl*tail_corr_ener;
11235 return scl*(tail_corr_virial + vdw_lambda_1*tail_corr_virial_1 +
11236 vdw_lambda_2*tail_corr_virial_2);
11239 return scl*tail_corr_virial;
11244 #ifdef MEM_OPT_VERSION 11247 int Molecule::checkExclByIdx(
int idx1,
int atom1,
int atom2)
const {
11249 int amin = exclChkSigPool[idx1].min;
11250 int amax = exclChkSigPool[idx1].max;
11251 int dist21 = atom2 - atom1;
11252 if ( dist21 < amin || dist21 > amax )
return 0;
11253 else return exclChkSigPool[idx1].flags[dist21-amin];
11259 int amin = all_exclusions[atom1].min;
11260 int amax = all_exclusions[atom1].max;
11261 if ( atom2 < amin || atom2 > amax )
return 0;
11262 else return all_exclusions[atom1].flags[atom2-amin];
11282 is_lonepairs_psf = 0;
11290 read_parm(amber_data);
11292 #ifndef MEM_OPT_VERSION 11295 assignLCPOTypes( 1 );
11307 is_lonepairs_psf = 0;
11315 read_parm(amber_data);
11317 #ifndef MEM_OPT_VERSION 11320 assignLCPOTypes( 1 );
11325 #ifdef MEM_OPT_VERSION 11326 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
11328 int i,j,ntheth,nphih,current_index,a1,a2,
11329 max,min,index,found;
11332 NAMD_die(
"No data read from parm file yet!");
11335 numAtoms = amber_data->
Natom;
11336 atoms =
new Atom[numAtoms];
11348 if (atoms == NULL || atomNames == NULL )
11349 NAMD_die(
"memory allocation failed when reading atom information");
11352 for (i = 0; i < numAtoms; ++i) {
11353 const int resname_length = amber_data->
ResNames[amber_data->
AtomRes[i]].size();
11354 const int atomname_length = amber_data->
AtomNames[i].size();
11355 const int atomtype_length = amber_data->
AtomSym[i].size();
11356 atomNames[i].resname = nameArena->getNewArray(resname_length+1);
11357 atomNames[i].atomname = nameArena->getNewArray(atomname_length+1);
11358 atomNames[i].atomtype = nameArena->getNewArray(atomtype_length+1);
11359 if (atomNames[i].resname == NULL || atomNames[i].atomname == NULL || atomNames[i].atomtype == NULL)
11360 NAMD_die(
"memory allocation failed when reading atom information");
11362 strncpy(atomNames[i].resname, amber_data->
ResNames[amber_data->
AtomRes[i]].c_str(), resname_length+1);
11363 strncpy(atomNames[i].atomname, amber_data->
AtomNames[i].c_str(), atomname_length+1);
11364 strncpy(atomNames[i].atomtype, amber_data->
AtomSym[i].c_str(), atomtype_length+1);
11366 strtok(atomNames[i].resname,
" ");
11367 strtok(atomNames[i].atomname,
" ");
11368 strtok(atomNames[i].atomtype,
" ");
11369 atoms[i].mass = amber_data->
Masses[i];
11416 atoms[i].charge = amber_data->
Charges[i] / 18.2223;
11417 atoms[i].vdw_type = amber_data->
Iac[i] - 1;
11419 if ( tmpResLookup ) tmpResLookup =
11422 if(atomSegResids) {
11424 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
11431 }
else if (atoms[i].mass <= 0.05) {
11433 ++numZeroMassAtoms;
11434 }
else if (atoms[i].mass < 1.0) {
11440 }
else if (atoms[i].mass <=3.5) {
11442 }
else if ((atomNames[i].atomname[0] ==
'O') &&
11443 (atoms[i].mass >= 14.0) &&
11444 (atoms[i].mass <= 18.0)) {
11457 if (amber_data->
Nbonh + amber_data->
Nbona > 0) {
11459 if (bonds == NULL || amber_data->
Nbona < 0)
11460 NAMD_die(
"memory allocation failed when reading bond information");
11462 for (i=0; i<amber_data->
Nbonh; ++i)
11463 { bonds[numBonds].atom1 = amber_data->
BondHAt1[i] / 3;
11464 bonds[numBonds].atom2 = amber_data->
BondHAt2[i] / 3;
11465 bonds[numBonds].bond_type = amber_data->
BondHNum[i] - 1;
11466 if (bonds[numBonds].atom1>=numAtoms || bonds[numBonds].atom2>=numAtoms ||
11467 bonds[numBonds].bond_type>=amber_data->
Numbnd)
11468 {
char err_msg[128];
11469 sprintf(err_msg,
"BOND (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11477 for (i=amber_data->
Nbonh; i<amber_data->Nbonh+amber_data->
Nbona; ++i)
11478 { bonds[numBonds].atom1 = amber_data->
BondAt1[i-amber_data->
Nbonh] / 3;
11479 bonds[numBonds].atom2 = amber_data->
BondAt2[i-amber_data->
Nbonh] / 3;
11480 bonds[numBonds].bond_type = amber_data->
BondNum[i-amber_data->
Nbonh] - 1;
11481 if (bonds[i].atom1>=numAtoms || bonds[i].atom2>=numAtoms ||
11482 bonds[i].bond_type>=amber_data->
Numbnd)
11483 {
char err_msg[128];
11484 sprintf(err_msg,
"BOND (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-amber_data->
Nbonh);
11493 if ( numBonds != amber_data->
Nbonh + amber_data->
Nbona) {
11495 " bonds with zero force constants.\n" <<
endi;
11497 "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
endi;
11503 { ntheth = amber_data->
Ntheth;
11504 angles =
new Angle[numAngles];
11505 if (angles == NULL || numAngles < ntheth)
11506 NAMD_die(
"memory allocation failed when reading angle information");
11508 for (i=0; i<ntheth; ++i)
11509 { angles[i].atom1 = amber_data->
AngleHAt1[i] / 3;
11510 angles[i].atom2 = amber_data->
AngleHAt2[i] / 3;
11511 angles[i].atom3 = amber_data->
AngleHAt3[i] / 3;
11512 angles[i].angle_type = amber_data->
AngleHNum[i] - 1;
11513 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11514 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11515 {
char err_msg[128];
11516 sprintf(err_msg,
"ANGLE (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11521 for (i=ntheth; i<numAngles; ++i)
11522 { angles[i].atom1 = amber_data->
AngleAt1[i-ntheth] / 3;
11523 angles[i].atom2 = amber_data->
AngleAt2[i-ntheth] / 3;
11524 angles[i].atom3 = amber_data->
AngleAt3[i-ntheth] / 3;
11525 angles[i].angle_type = amber_data->
AngleNum[i-ntheth] - 1;
11526 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11527 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11528 {
char err_msg[128];
11529 sprintf(err_msg,
"ANGLE (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-ntheth);
11544 if (exclusions == NULL && amber_data->
Nnb > 0)
11545 NAMD_die(
"memory allocation failed when reading exclusion information");
11547 for (i=0; i<numAtoms; ++i)
11548 for (j=0; j<amber_data->
Iblo[i]; ++j)
11549 {
if (current_index >= amber_data->
Nnb)
11550 {
char err_msg[128];
11551 sprintf(err_msg,
"EXCLUSION INDEX EXCEEDS NUMBER OF EXLCUSIONS %d IN AMBER FILE, AT ATOM #%d\n",
11552 amber_data->
Nnb, i+1);
11557 if (amber_data->
ExclAt[current_index] != 0)
11560 a2 = amber_data->
ExclAt[current_index] - 1;
11566 sprintf(err_msg,
"Atom #%d has exclusion with atom #%d, in reverse order.", i+1, a2+1);
11570 {
char err_msg[128];
11571 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN AMBER FILE\n", i+1);
11574 else if (a2 >= numAtoms)
11575 {
char err_msg[128];
11576 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN AMBER FILE",
11577 a2+1, numAtoms, current_index+1);
11580 exclusions[numExclusions].atom1 = i;
11581 exclusions[numExclusions].atom2 = a2;
11586 if (current_index < amber_data->Nnb)
11587 {
char err_msg[128];
11588 sprintf(err_msg,
"Num of exclusions recorded (%d) is smaller than what it's supposed to be (%d)",
11589 current_index,amber_data->
Nnb);
11595 numDihedrals = amber_data->
Nphih + amber_data->
Nphia;
11596 if (numDihedrals > 0)
11597 { nphih = amber_data->
Nphih;
11598 dihedrals =
new Dihedral[numDihedrals];
11599 if (dihedrals == NULL || numDihedrals < nphih)
11600 NAMD_die(
"memory allocation failed when reading dihedral information");
11602 for (i=0; i<nphih; ++i)
11603 { dihedrals[i].atom1 = amber_data->
DihHAt1[i] / 3;
11604 dihedrals[i].atom2 = amber_data->
DihHAt2[i] / 3;
11605 dihedrals[i].atom3 = amber_data->
DihHAt3[i] / 3;
11606 dihedrals[i].atom4 = amber_data->
DihHAt4[i] / 3;
11607 dihedrals[i].dihedral_type = amber_data->
DihHNum[i] - 1;
11610 for (i=nphih; i<numDihedrals; ++i)
11611 { dihedrals[i].atom1 = amber_data->
DihAt1[i-nphih] / 3;
11612 dihedrals[i].atom2 = amber_data->
DihAt2[i-nphih] / 3;
11613 dihedrals[i].atom3 = amber_data->
DihAt3[i-nphih] / 3;
11614 dihedrals[i].atom4 = amber_data->
DihAt4[i-nphih] / 3;
11615 dihedrals[i].dihedral_type = amber_data->
DihNum[i-nphih] - 1;
11630 for (i=0; i<numDihedrals; ++i)
11631 {
if (dihedrals[i].atom3 < 0 || dihedrals[i].atom4 < 0)
11632 { dihedrals[i].atom3 = abs(dihedrals[i].atom3);
11633 dihedrals[i].atom4 = abs(dihedrals[i].atom4);
11636 {
if (dihedrals[i].atom1 < dihedrals[i].atom4)
11637 a1=dihedrals[i].atom1, a2=dihedrals[i].atom4;
11639 a1=dihedrals[i].atom4, a2=dihedrals[i].atom1;
11643 min=0, max=numExclusions-1;
11644 while (!found && min<=max)
11645 { index = (min+max)/2;
11646 if (exclusions[index].atom1 == a1)
11648 else if (exclusions[index].atom1 < a1)
11654 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11657 for (j=index-1; j>=0 && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; --j);
11658 if (j<0 || exclusions[j].atom1!=a1)
11659 for (j=index; j<numExclusions && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; ++j);
11660 if (j<numExclusions && exclusions[j].atom1==a1)
11661 exclusions[j].modified = 1;
11663 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11665 if (dihedrals[i].atom1>=numAtoms || dihedrals[i].atom2>=numAtoms ||
11666 dihedrals[i].atom3>=numAtoms || dihedrals[i].atom4>=numAtoms ||
11667 dihedrals[i].dihedral_type>=amber_data->
Nptra)
11668 {
char err_msg[128];
11669 sprintf(err_msg,
"DIHEDRAL # %d OVERFLOW IN PARM FILE", i+1);
11675 if (amber_data -> HasCMAP) {
11677 if (numCrossterms > 0) {
11678 crossterms =
new Crossterm[numCrossterms];
11680 for (i = 0; i < numCrossterms; ++i) {
11683 crossterms[i].atom2 = amber_data->
CMAPIndex[6*i+1]-1;
11684 crossterms[i].atom3 = amber_data->
CMAPIndex[6*i+2]-1;
11685 crossterms[i].atom4 = amber_data->
CMAPIndex[6*i+3]-1;
11687 crossterms[i].atom5 = amber_data->
CMAPIndex[6*i+1]-1;
11688 crossterms[i].atom6 = amber_data->
CMAPIndex[6*i+2]-1;
11689 crossterms[i].atom7 = amber_data->
CMAPIndex[6*i+3]-1;
11690 crossterms[i].atom8 = amber_data->
CMAPIndex[6*i+4]-1;
11692 crossterms[i].crossterm_type = amber_data->
CMAPIndex[6*i+5]-1;
11697 numRealBonds = numBonds;
11698 build_atom_status();
11714 void Molecule::read_parm(
Ambertoppar *amber_data)
11716 #ifdef MEM_OPT_VERSION 11717 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
11719 int i,j,ntheth,nphih,current_index,a1,a2,
11720 max,min,index,found;
11723 NAMD_die(
"No data read from parm file yet!");
11726 numAtoms = amber_data->
Natom;
11727 atoms =
new Atom[numAtoms];
11734 if (atoms == NULL || atomNames == NULL )
11735 NAMD_die(
"memory allocation failed when reading atom information");
11737 for (i=0; i<numAtoms; ++i)
11738 { atomNames[i].resname = nameArena->getNewArray(5);
11739 atomNames[i].atomname = nameArena->getNewArray(5);
11740 atomNames[i].atomtype = nameArena->getNewArray(5);
11741 if (atomNames[i].resname == NULL || atomNames[i].atomname == NULL || atomNames[i].atomtype == NULL)
11742 NAMD_die(
"memory allocation failed when reading atom information");
11743 for (j=0; j<4; ++j)
11744 { atomNames[i].resname[j] = amber_data->
ResNames[amber_data->
AtomRes[i]*4+j];
11745 atomNames[i].atomname[j] = amber_data->
AtomNames[i*4+j];
11746 atomNames[i].atomtype[j] = amber_data->
AtomSym[i*4+j];
11748 atomNames[i].resname[4] = atomNames[i].atomname[4] = atomNames[i].atomtype[4] =
'\0';
11749 strtok(atomNames[i].resname,
" ");
11750 strtok(atomNames[i].atomname,
" ");
11751 strtok(atomNames[i].atomtype,
" ");
11752 atoms[i].mass = amber_data->
Masses[i];
11754 atoms[i].charge = amber_data->
Charges[i] / 18.2223;
11755 atoms[i].vdw_type = amber_data->
Iac[i] - 1;
11758 if ( tmpResLookup ) tmpResLookup =
11761 if(atomSegResids) {
11763 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
11771 }
else if (atoms[i].mass <= 0.05) {
11772 ++numZeroMassAtoms;
11774 }
else if (atoms[i].mass < 1.0) {
11776 }
else if (atoms[i].mass <=3.5) {
11778 }
else if ((atomNames[i].atomname[0] ==
'O') &&
11779 (atoms[i].mass >= 14.0) &&
11780 (atoms[i].mass <= 18.0)) {
11793 if (amber_data->
Nbonh + amber_data->
Nbona > 0)
11795 if (bonds == NULL || amber_data->
Nbona < 0)
11796 NAMD_die(
"memory allocation failed when reading bond information");
11798 for (i=0; i<amber_data->
Nbonh; ++i)
11799 { bonds[numBonds].atom1 = amber_data->
BondHAt1[i] / 3;
11800 bonds[numBonds].atom2 = amber_data->
BondHAt2[i] / 3;
11801 bonds[numBonds].bond_type = amber_data->
BondHNum[i] - 1;
11802 if (bonds[numBonds].atom1>=numAtoms || bonds[numBonds].atom2>=numAtoms ||
11803 bonds[numBonds].bond_type>=amber_data->
Numbnd)
11804 {
char err_msg[128];
11805 sprintf(err_msg,
"BOND (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11813 for (i=amber_data->
Nbonh; i<amber_data->Nbonh+amber_data->
Nbona; ++i)
11814 { bonds[numBonds].atom1 = amber_data->
BondAt1[i-amber_data->
Nbonh] / 3;
11815 bonds[numBonds].atom2 = amber_data->
BondAt2[i-amber_data->
Nbonh] / 3;
11816 bonds[numBonds].bond_type = amber_data->
BondNum[i-amber_data->
Nbonh] - 1;
11817 if (bonds[i].atom1>=numAtoms || bonds[i].atom2>=numAtoms ||
11818 bonds[i].bond_type>=amber_data->
Numbnd)
11819 {
char err_msg[128];
11820 sprintf(err_msg,
"BOND (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-amber_data->
Nbonh);
11829 if ( numBonds != amber_data->
Nbonh + amber_data->
Nbona) {
11831 " bonds with zero force constants.\n" <<
endi;
11833 "Will get H-H distance in rigid H2O from H-O-H angle.\n" <<
endi;
11839 { ntheth = amber_data->
Ntheth;
11840 angles =
new Angle[numAngles];
11841 if (angles == NULL || numAngles < ntheth)
11842 NAMD_die(
"memory allocation failed when reading angle information");
11844 for (i=0; i<ntheth; ++i)
11845 { angles[i].atom1 = amber_data->
AngleHAt1[i] / 3;
11846 angles[i].atom2 = amber_data->
AngleHAt2[i] / 3;
11847 angles[i].atom3 = amber_data->
AngleHAt3[i] / 3;
11848 angles[i].angle_type = amber_data->
AngleHNum[i] - 1;
11849 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11850 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11851 {
char err_msg[128];
11852 sprintf(err_msg,
"ANGLE (WITH H) # %d OVERFLOW IN PARM FILE", i+1);
11857 for (i=ntheth; i<numAngles; ++i)
11858 { angles[i].atom1 = amber_data->
AngleAt1[i-ntheth] / 3;
11859 angles[i].atom2 = amber_data->
AngleAt2[i-ntheth] / 3;
11860 angles[i].atom3 = amber_data->
AngleAt3[i-ntheth] / 3;
11861 angles[i].angle_type = amber_data->
AngleNum[i-ntheth] - 1;
11862 if (angles[i].atom1>=numAtoms || angles[i].atom2>=numAtoms ||
11863 angles[i].atom3>=numAtoms || angles[i].angle_type>=amber_data->
Numang)
11864 {
char err_msg[128];
11865 sprintf(err_msg,
"ANGLE (WITHOUT H) # %d OVERFLOW IN PARM FILE", i+1-ntheth);
11880 if (exclusions == NULL && amber_data->
Nnb > 0)
11881 NAMD_die(
"memory allocation failed when reading exclusion information");
11883 for (i=0; i<numAtoms; ++i)
11884 for (j=0; j<amber_data->
Iblo[i]; ++j)
11885 {
if (current_index >= amber_data->
Nnb)
11886 {
char err_msg[128];
11887 sprintf(err_msg,
"EXCLUSION INDEX EXCEEDS NUMBER OF EXLCUSIONS %d IN AMBER FILE, AT ATOM #%d\n",
11888 amber_data->
Nnb, i+1);
11893 if (amber_data->
ExclAt[current_index] != 0)
11896 a2 = amber_data->
ExclAt[current_index] - 1;
11902 sprintf(err_msg,
"Atom #%d has exclusion with atom #%d, in reverse order.", i+1, a2+1);
11906 {
char err_msg[128];
11907 sprintf(err_msg,
"ATOM %d EXCLUDED FROM ITSELF IN AMBER FILE\n", i+1);
11910 else if (a2 >= numAtoms)
11911 {
char err_msg[128];
11912 sprintf(err_msg,
"EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN AMBER FILE",
11913 a2+1, numAtoms, current_index+1);
11916 exclusions[numExclusions].atom1 = i;
11917 exclusions[numExclusions].atom2 = a2;
11922 if (current_index < amber_data->Nnb)
11923 {
char err_msg[128];
11924 sprintf(err_msg,
"Num of exclusions recorded (%d) is smaller than what it's supposed to be (%d)",
11925 current_index,amber_data->
Nnb);
11931 numDihedrals = amber_data->
Nphih + amber_data->
Nphia;
11932 if (numDihedrals > 0)
11933 { nphih = amber_data->
Nphih;
11934 dihedrals =
new Dihedral[numDihedrals];
11935 if (dihedrals == NULL || numDihedrals < nphih)
11936 NAMD_die(
"memory allocation failed when reading dihedral information");
11938 for (i=0; i<nphih; ++i)
11939 { dihedrals[i].atom1 = amber_data->
DihHAt1[i] / 3;
11940 dihedrals[i].atom2 = amber_data->
DihHAt2[i] / 3;
11941 dihedrals[i].atom3 = amber_data->
DihHAt3[i] / 3;
11942 dihedrals[i].atom4 = amber_data->
DihHAt4[i] / 3;
11943 dihedrals[i].dihedral_type = amber_data->
DihHNum[i] - 1;
11946 for (i=nphih; i<numDihedrals; ++i)
11947 { dihedrals[i].atom1 = amber_data->
DihAt1[i-nphih] / 3;
11948 dihedrals[i].atom2 = amber_data->
DihAt2[i-nphih] / 3;
11949 dihedrals[i].atom3 = amber_data->
DihAt3[i-nphih] / 3;
11950 dihedrals[i].atom4 = amber_data->
DihAt4[i-nphih] / 3;
11951 dihedrals[i].dihedral_type = amber_data->
DihNum[i-nphih] - 1;
11966 for (i=0; i<numDihedrals; ++i)
11967 {
if (dihedrals[i].atom3 < 0 || dihedrals[i].atom4 < 0)
11968 { dihedrals[i].atom3 = abs(dihedrals[i].atom3);
11969 dihedrals[i].atom4 = abs(dihedrals[i].atom4);
11972 {
if (dihedrals[i].atom1 < dihedrals[i].atom4)
11973 a1=dihedrals[i].atom1, a2=dihedrals[i].atom4;
11975 a1=dihedrals[i].atom4, a2=dihedrals[i].atom1;
11979 min=0, max=numExclusions-1;
11980 while (!found && min<=max)
11981 { index = (min+max)/2;
11982 if (exclusions[index].atom1 == a1)
11984 else if (exclusions[index].atom1 < a1)
11990 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
11993 for (j=index-1; j>=0 && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; --j);
11994 if (j<0 || exclusions[j].atom1!=a1)
11995 for (j=index; j<numExclusions && exclusions[j].atom2!=a2 && exclusions[j].atom1==a1; ++j);
11996 if (j<numExclusions && exclusions[j].atom1==a1)
11997 exclusions[j].modified = 1;
11999 NAMD_die(
"1-4 interaction in dihedral not found in exclusion list!");
12001 if (dihedrals[i].atom1>=numAtoms || dihedrals[i].atom2>=numAtoms ||
12002 dihedrals[i].atom3>=numAtoms || dihedrals[i].atom4>=numAtoms ||
12003 dihedrals[i].dihedral_type>=amber_data->
Nptra)
12004 {
char err_msg[128];
12005 sprintf(err_msg,
"DIHEDRAL # %d OVERFLOW IN PARM FILE", i+1);
12011 numRealBonds = numBonds;
12012 build_atom_status();
12031 read_parm(gromacsTopFile);
12033 #ifndef MEM_OPT_VERSION 12036 assignLCPOTypes( 3 );
12054 #ifdef MEM_OPT_VERSION 12055 NAMD_die(
"When reading a compressed file or using the memory-optimized version, amber data is not supported!");
12063 atoms =
new Atom[numAtoms];
12070 if (atoms == NULL || atomNames == NULL )
12071 NAMD_die(
"memory allocation failed when reading atom information");
12075 for (i=0; i<numAtoms; ++i) {
12076 char *resname = nameArena->getNewArray(11);
12077 char *atomname = nameArena->getNewArray(11);
12078 char *atomtype = nameArena->getNewArray(11);
12079 int resnum,typenum;
12082 if (resname == NULL || atomname == NULL || atomtype == NULL)
12083 NAMD_die(
"memory allocation failed when reading atom information");
12086 gf->
getAtom(i,&resnum,resname,atomname,atomtype,&typenum,
12089 atomNames[i].resname = resname;
12090 atomNames[i].atomname = atomname;
12091 atomNames[i].atomtype = atomtype;
12092 atoms[i].mass = mass;
12093 atoms[i].charge = charge;
12094 atoms[i].vdw_type = typenum;
12097 if ( tmpResLookup ) tmpResLookup =
12098 tmpResLookup->
append(
"MAIN", resnum+1, i);
12100 if(atomSegResids) {
12102 memcpy(one->
segname,
"MAIN", strlen(
"MAIN")+1);
12103 one->
resid = resnum+1;
12112 }
else if (atoms[i].mass <= 0.05) {
12114 }
else if (atoms[i].mass < 1.0) {
12116 }
else if (atoms[i].mass <=3.5) {
12118 }
else if ((atomNames[i].atomname[0] ==
'O') &&
12119 (atoms[i].mass >= 14.0) &&
12120 (atoms[i].mass <= 18.0)) {
12127 bonds =
new Bond[numBonds];
12129 NAMD_die(
"memory allocation failed when reading bond information");
12130 for(i=0;i<numBonds;i++) {
12133 gf->
getBond(i,&atom1,&atom2,&type);
12134 bonds[i].atom1 = atom1;
12135 bonds[i].atom2 = atom2;
12136 bonds[i].bond_type = (
Index)type;
12141 angles =
new Angle[numAngles];
12142 if (angles == NULL)
12143 NAMD_die(
"memory allocation failed when reading angle information");
12144 for(i=0;i<numAngles;i++) {
12146 int atom1,atom2,atom3;
12147 gf->
getAngle(i,&atom1,&atom2,&atom3,&type);
12149 angles[i].atom1 = atom1;
12150 angles[i].atom2 = atom2;
12151 angles[i].atom3 = atom3;
12153 angles[i].angle_type=type;
12157 exclusions =
new Exclusion[numExclusions];
12221 dihedrals =
new Dihedral[numDihedrals];
12222 if (dihedrals == NULL)
12223 NAMD_die(
"memory allocation failed when reading dihedral information");
12224 for(i=0;i<numDihedrals;i++) {
12226 int atom1,atom2,atom3,atom4;
12227 gf->
getDihedral(i,&atom1,&atom2,&atom3,&atom4,&type);
12228 dihedrals[i].atom1 = atom1;
12229 dihedrals[i].atom2 = atom2;
12230 dihedrals[i].atom3 = atom3;
12231 dihedrals[i].atom4 = atom4;
12232 dihedrals[i].dihedral_type = type;
12245 const_cast<GromacsTopFile*
>(gf)->getPairLJArrays2(indxLJA, indxLJB, pairC6, pairC12);
12248 gromacsPair_type[i] = i;
12257 pointerToLJBeg =
new int[numAtoms];
12258 pointerToLJEnd =
new int[numAtoms];
12260 for(
int i=0; i < numAtoms; i++) {
12261 pointerToLJBeg[i] = -1;
12262 pointerToLJEnd[i] = -2;
12265 if(pointerToLJBeg[indxLJA[i]] == -1) {
12266 pointerToLJBeg[indxLJA[i]] = i;
12267 oldIndex = indxLJA[i];
12269 pointerToLJEnd[oldIndex] = i;
12282 const_cast<GromacsTopFile*
>(gf)->getPairGaussArrays2(indxGaussA, indxGaussB, gA, gMu1, giSigma1, gMu2, giSigma2, gRepulsive);
12285 pointerToGaussBeg =
new int[numAtoms];
12286 pointerToGaussEnd =
new int[numAtoms];
12287 for(
int i=0; i < numAtoms; i++) {
12288 pointerToGaussBeg[i] = -1;
12289 pointerToGaussEnd[i] = -2;
12293 if(pointerToGaussBeg[indxGaussA[i]] == -1) {
12294 pointerToGaussBeg[indxGaussA[i]] = i;
12295 oldIndex = indxGaussA[i];
12297 pointerToGaussEnd[oldIndex] = i;
12300 iout <<
iINFO <<
"Finished reading explicit pair from Gromacs file:\n" <<
12301 iINFO <<
"Found a total of: " << numPair <<
" explicit pairs--of which: " <<
12303 " are Gaussian style pairs.\n" <<
endi;
12307 #if GROMACS_EXCLUSIONS 12310 int* atom1 =
new int[numExclusions];
12311 int* atom2 =
new int[numExclusions];
12312 for(
int j=0; j<numExclusions;j++) {
12318 read_exclusions(atom1,atom2,numExclusions);
12382 numRealBonds = numBonds;
12383 build_atom_status();
12388 #ifndef MEM_OPT_VERSION 12402 #ifdef MEM_OPT_VERSION 12404 Index Molecule::insert_new_mass(
Real newMass){
12406 for(
int i=massPoolSize-1; i>=0; i--){
12407 if(fabs(atomMassPool[i]-newMass)<=1e-6)
12411 Real *tmp =
new Real[massPoolSize+1];
12412 tmp[massPoolSize] = newMass;
12413 memcpy((
void *)tmp, (
const void *)atomMassPool,
sizeof(
Real)*massPoolSize);
12414 delete [] atomMassPool;
12415 atomMassPool = tmp;
12417 return (
Index)(massPoolSize-1);
12420 void Molecule::addNewExclSigPool(
const vector<ExclusionSignature>& newExclSigPool){
12422 for(
int i=0; i<exclSigPoolSize; i++)
12423 tmpExclSigPool[i] = exclSigPool[i];
12424 for(
int i=0; i<newExclSigPool.size(); i++)
12425 tmpExclSigPool[i+exclSigPoolSize] = newExclSigPool[i];
12427 exclSigPoolSize += newExclSigPool.size();
12428 exclSigPool = tmpExclSigPool;
12432 msg->
put((
short)tupleType);
12433 msg->
put(numOffset);
12434 msg->
put(numOffset, offset);
12435 msg->
put(tupleParamType);
12444 msg->
get(numOffset);
12446 offset =
new int[numOffset];
12447 msg->
get(numOffset*
sizeof(
int), (
char *)offset);
12449 msg->
get(tupleParamType);
12455 for(
int i=0; i<bondCnt; i++)
12456 bondSigs[i].pack(msg);
12458 msg->
put(angleCnt);
12459 for(
int i=0; i<angleCnt; i++)
12460 angleSigs[i].pack(msg);
12462 msg->
put(dihedralCnt);
12463 for(
int i=0; i<dihedralCnt; i++)
12464 dihedralSigs[i].pack(msg);
12466 msg->
put(improperCnt);
12467 for(
int i=0; i<improperCnt; i++)
12468 improperSigs[i].pack(msg);
12470 msg->
put(crosstermCnt);
12471 for(
int i=0; i<crosstermCnt; i++)
12472 crosstermSigs[i].pack(msg);
12475 msg->
put(gromacsPairCnt);
12476 for(
int i=0; i<gromacsPairCnt; i++)
12477 gromacsPairSigs[i].pack(msg);
12482 delete [] bondSigs;
12485 for(
int i=0; i<bondCnt; i++)
12486 bondSigs[i].unpack(msg);
12487 }
else bondSigs = NULL;
12489 msg->
get(angleCnt);
12490 delete [] angleSigs;
12493 for(
int i=0; i<angleCnt; i++)
12494 angleSigs[i].unpack(msg);
12495 }
else angleSigs = NULL;
12497 msg->
get(dihedralCnt);
12498 delete [] dihedralSigs;
12501 for(
int i=0; i<dihedralCnt; i++)
12502 dihedralSigs[i].unpack(msg);
12503 }
else dihedralSigs = NULL;
12505 msg->
get(improperCnt);
12506 delete [] improperSigs;
12509 for(
int i=0; i<improperCnt; i++)
12510 improperSigs[i].unpack(msg);
12511 }
else improperSigs = NULL;
12513 msg->
get(crosstermCnt);
12514 delete [] crosstermSigs;
12515 if(crosstermCnt>0){
12517 for(
int i=0; i<crosstermCnt; i++)
12518 crosstermSigs[i].unpack(msg);
12519 }
else crosstermSigs = NULL;
12523 msg->
get(gromacsPairCnt);
12524 delete [] gromacsPairSigs;
12525 if(gromacsPairCnt>0){
12527 for(
int i=0; i<gromacsPairCnt; i++)
12528 gromacsPairSigs[i].unpack(msg);
12529 }
else gromacsPairSigs = NULL;
12543 origTupleCnt = bondCnt;
12544 tupleSigs= bondSigs;
12545 for(
int i=0; i<origTupleCnt; i++){
12550 delete [] tupleSigs;
12552 }
else if(bondCnt!=origTupleCnt){
12555 for(
int i=0; i<origTupleCnt; i++){
12557 newTupleSigs[idx] = tupleSigs[i];
12561 delete [] tupleSigs;
12562 bondSigs = newTupleSigs;
12568 origTupleCnt = angleCnt;
12569 tupleSigs = angleSigs;
12570 for(
int i=0; i<origTupleCnt; i++){
12575 delete [] tupleSigs;
12577 }
else if(angleCnt!=origTupleCnt){
12580 for(
int i=0; i<origTupleCnt; i++){
12582 newTupleSigs[idx] = tupleSigs[i];
12586 delete [] tupleSigs;
12587 angleSigs = newTupleSigs;
12593 origTupleCnt = dihedralCnt;
12594 tupleSigs = dihedralSigs;
12595 for(
int i=0; i<origTupleCnt; i++){
12599 if(dihedralCnt==0){
12600 delete [] tupleSigs;
12601 dihedralSigs = NULL;
12602 }
else if(dihedralCnt!=origTupleCnt){
12605 for(
int i=0; i<origTupleCnt; i++){
12607 newTupleSigs[idx] = tupleSigs[i];
12611 delete [] tupleSigs;
12612 dihedralSigs = newTupleSigs;
12619 origTupleCnt = improperCnt;
12620 tupleSigs = improperSigs;
12621 for(
int i=0; i<origTupleCnt; i++){
12625 if(improperCnt==0){
12626 delete [] tupleSigs;
12627 improperSigs = NULL;
12628 }
else if(improperCnt!=origTupleCnt){
12631 for(
int i=0; i<origTupleCnt; i++){
12633 newTupleSigs[idx] = tupleSigs[i];
12637 delete [] tupleSigs;
12638 improperSigs = newTupleSigs;
12644 origTupleCnt = crosstermCnt;
12645 tupleSigs = crosstermSigs;
12646 for(
int i=0; i<origTupleCnt; i++){
12650 if(crosstermCnt==0){
12651 delete [] tupleSigs;
12652 crosstermSigs = NULL;
12653 }
else if(crosstermCnt!=origTupleCnt){
12656 for(
int i=0; i<origTupleCnt; i++){
12658 newTupleSigs[idx] = tupleSigs[i];
12662 delete [] tupleSigs;
12663 crosstermSigs = newTupleSigs;
12670 origTupleCnt = gromacsPairCnt;
12671 tupleSigs = gromacsPairSigs;
12672 for(
int i=0; i<origTupleCnt; i++){
12676 if(gromacsPairCnt==0){
12677 delete [] tupleSigs;
12678 gromacsPairSigs = NULL;
12679 }
else if(gromacsPairCnt!=origTupleCnt){
12682 for(
int i=0; i<origTupleCnt; i++){
12684 newTupleSigs[idx] = tupleSigs[i];
12688 delete [] tupleSigs;
12689 gromacsPairSigs = newTupleSigs;
12699 for(
int i=0; i<fullExclCnt; i++){
12700 if(fullOffset[i]==0)
continue;
12705 delete [] fullOffset;
12707 }
else if(newCnt!=fullExclCnt){
12708 int *tmpOffset =
new int[newCnt];
12710 for(
int i=0; i<fullExclCnt; i++){
12711 if(fullOffset[i]==0)
continue;
12712 tmpOffset[newCnt] = fullOffset[i];
12715 delete [] fullOffset;
12716 fullOffset = tmpOffset;
12717 fullExclCnt = newCnt;
12722 for(
int i=0; i<modExclCnt; i++){
12723 if(modOffset[i]==0)
continue;
12728 delete [] modOffset;
12730 }
else if(newCnt!=modExclCnt){
12731 int *tmpOffset =
new int[newCnt];
12733 for(
int i=0; i<modExclCnt; i++){
12734 if(modOffset[i]==0)
continue;
12735 tmpOffset[newCnt] = modOffset[i];
12738 delete [] modOffset;
12739 modOffset = tmpOffset;
12740 modExclCnt = newCnt;
12755 int high = fullExclCnt-1;
12756 int mid = (low+high)/2;
12758 if(offset<fullOffset[mid]){
12760 mid = (high+low)/2;
12761 }
else if(offset>fullOffset[mid]){
12763 mid = (high+low)/2;
12769 if(retidx!=-1)
return retidx;
12773 high = modExclCnt-1;
12774 mid = (low+high)/2;
12776 if(offset<modOffset[mid]){
12778 mid = (high+low)/2;
12779 }
else if(offset>modOffset[mid]){
12781 mid = (high+low)/2;
12791 msg->
put(fullExclCnt);
12792 msg->
put(fullExclCnt, fullOffset);
12793 msg->
put(modExclCnt);
12794 msg->
put(modExclCnt, modOffset);
12798 msg->
get(fullExclCnt);
12799 delete [] fullOffset;
12800 fullOffset =
new int[fullExclCnt];
12801 msg->
get(fullExclCnt*
sizeof(
int), (
char *)fullOffset);
12802 msg->
get(modExclCnt);
12803 delete [] modOffset;
12804 modOffset =
new int[modExclCnt];
12805 msg->
get(modExclCnt*
sizeof(
int), (
char *)modOffset);
12806 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 12812 #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)
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
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 *)