30 #define MIN_DEBUG_LEVEL 3
38 #ifndef CODE_REDUNDANT
39 #define CODE_REDUNDANT 0
60 const char* newSegName,
Real newQMID) {
65 atmIDs.push_back(newBegAtm) ;
73 std::vector<std::string>
split(
const std::string &text, std::string delimiter) {
75 std::vector<std::string> tokens;
76 std::size_t start = 0, end = 0;
78 while ((end = text.find(delimiter, start)) != std::string::npos) {
80 std::string temp = text.substr(start, end - start);
82 if (! temp.empty()) tokens.push_back(temp);
84 start = end + delimiter.length();
88 std::string temp = text.substr(start);
89 if (! temp.empty()) tokens.push_back(temp);
111 #ifdef MEM_OPT_VERSION
112 NAMD_die(
"QMMM interface is not supported in memory optimized builds.");
121 iout <<
iINFO <<
"Using the following PDB file for QM parameters: " <<
122 pdbFileName <<
"\n" <<
endi;
124 pdbP =
new PDB(pdbFileName);
126 iout <<
"pdbFile name not defined!\n\n" <<
endi ;
130 NAMD_die(
"Number of atoms in QM parameter PDB file doesn't match coordinate PDB");
133 qmElementArray =
new char*[
numAtoms] ;
140 std::set<Real> atmGrpSet;
142 std::vector<Real> grpChrgVec;
145 std::vector<BigReal> qmBondValVec;
147 std::vector<Real> qmGroupIDsVec;
150 std::map<Real,int> qmGrpIDMap ;
152 std::vector<int> qmAtmIndxVec, qmGrpSizeVec;
155 std::vector<std::string> strVec;
157 qmNoPC = simParams->
qmNoPC;
169 qmLSSTotalNumAtms = 0;
171 std::map<Real,std::vector<int> > lssGrpRefIDs;
173 int totalNumRefAtms = 0;
174 int lssClassicResIndx = 0 ;
175 int lssCurrClassResID = -1 ;
176 char lssCurrClassResSegID[5];
178 DebugM(4,
"LSS is ON! Processing QM solvent.\n") ;
180 if (resLookup == NULL)
181 NAMD_die(
"We need residue data to conduct LSS.");
186 for ( ; current; current = current->
next ) {
188 strVec =
split( std::string(current->
data) ,
" ");
190 if (strVec.size() != 3 ) {
191 iout <<
iERROR <<
"Format error in QM LSS size: "
194 NAMD_die(
"Error processing QM information.");
197 std::stringstream storConv ;
199 storConv << strVec[0] ;
202 if (storConv.fail()) {
203 iout <<
iERROR <<
"Error parsing QMLSSRef selection: "
206 NAMD_die(
"Error processing QM information.");
209 std::string segName = strVec[1].substr(0,4);
212 storConv << strVec[2];
215 if (storConv.fail()) {
216 iout <<
iERROR <<
"Error parsing QMLSSRef selection: "
219 NAMD_die(
"Error processing QM information.");
222 auto it = lssRefUsrSel.find(grpID) ;
223 if (it == lssRefUsrSel.end())
226 lssRefUsrSel[grpID].push_back(
refSelStr(segName, resID));
229 for (
auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
230 iout <<
iINFO <<
"LSS user selection for COM composition of group "
231 << it->first <<
":\n" <<
endi ;
232 for (
int i=0; i<it->second.size();i++) {
233 iout <<
iINFO <<
"Segment " << it->second[i].segid
234 <<
" ; residue " << it->second[i].resid <<
"\n" <<
endi ;
246 for (
int atmInd = 0 ; atmInd <
numAtoms; atmInd++) {
252 if ( strcmp(simParams->
qmColumn,
"beta") == 0 ){
263 qmBondValVec.push_back(bondVal);
267 if ( strcmp(simParams->
qmColumn,
"beta") == 0 ){
277 qmAtomGroup[atmInd] = atmGrp;
281 qmElementArray[atmInd] =
new char[3];
282 strncpy(qmElementArray[atmInd],pdbP->
atom(atmInd)->
element(),3);
289 if ( fixedAtomFlags[atmInd] == 1 ) {
290 iout <<
iERROR <<
"QM atom cannot be fixed in space!\n"
292 NAMD_die(
"Error processing QM information.");
299 if (simParams->
qmVDW) {
301 std::string qmType(qmElementArray[atmInd]) ;
304 std::remove_if(qmType.begin(), qmType.end(), isspace ),
307 qmType = std::string(
"q") + qmType;
324 qmAtmIndxVec.push_back(atmInd);
326 auto retVal = atmGrpSet.insert(atmGrp);
333 qmGrpIDMap.insert(std::pair<BigReal,int>(atmGrp, qmNumGrps));
336 qmGroupIDsVec.push_back(atmGrp);
344 qmGrpSizeVec.push_back(1);
348 grpChrgVec.push_back( atoms[atmInd].
charge );
353 qmGrpSizeVec[ qmGrpIDMap[atmGrp] ]++ ;
356 grpChrgVec[ qmGrpIDMap[atmGrp] ] += atoms[atmInd].
charge;
370 resLookup->
lookup(segName,resID,&i, &end);
379 resP->
atmIDs.push_back(atmInd) ;
398 auto itNewGrp = lssGrpRefIDs.find(atmGrp) ;
399 if (itNewGrp == lssGrpRefIDs.end()) {
400 lssGrpRefIDs.insert(std::pair<
Real, std::vector<int> >(
401 atmGrp, std::vector<int>() ) );
411 for(
int i=0; i<lssRefUsrSel[atmGrp].size(); i++) {
412 if (lssRefUsrSel[atmGrp][i].resid == resID &&
413 strncmp(lssRefUsrSel[atmGrp][i].segid.c_str(),
416 lssGrpRefIDs[atmGrp].push_back(atmInd);
426 lssGrpRefIDs[atmGrp].push_back(atmInd);
437 else if (atmGrp == 0) {
448 if (lssCurrClassResID < 0) {
449 lssCurrClassResID = resID ;
450 strncpy(lssCurrClassResSegID, segName,5);
451 lssClassicResIndx = 0;
453 else if (lssCurrClassResID != resID ||
454 strcmp(lssCurrClassResSegID, segName) != 0 ) {
455 lssCurrClassResID = resID ;
456 strncpy(lssCurrClassResSegID, segName,5);
460 qmClassicSolv.insert(std::pair<int,int>(atmInd,lssClassicResIndx));
465 else if(atmGrp < 0) {
466 iout <<
iERROR <<
"QM group ID cannot be less than zero!\n" <<
endi;
467 NAMD_die(
"Error processing QM information.");
473 if (lssGrpRes.
size() == 0)
474 NAMD_die(
"LSS was activated but there are no QM solvent molecules!\n");
476 for (
auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
477 auto itTarget = qmGrpIDMap.find(it->first);
478 if (itTarget == qmGrpIDMap.end()) {
479 iout <<
iERROR <<
"QM group ID for LSS could not be found in input!"
480 <<
" QM group ID: " << it->first <<
"\n" <<
endi;
481 NAMD_die(
"Error processing QM information.");
485 DebugM(3,
"We found " << lssClassicResIndx <<
" classical solvent molecules totalling "
486 << qmClassicSolv.size() <<
" atoms.\n" );
490 qmNumQMAtoms = qmAtmIndxVec.size();
492 if (qmNumQMAtoms == 0)
493 NAMD_die(
"No QM atoms were found in this QM simulation!") ;
495 iout <<
iINFO <<
"Number of QM atoms (excluding Dummy atoms): " <<
496 qmNumQMAtoms <<
"\n" <<
endi;
498 qmAtmChrg =
new Real[qmNumQMAtoms] ;
499 qmAtmIndx =
new int[qmNumQMAtoms] ;
500 for (
int i=0; i<qmNumQMAtoms; i++) {
504 qmAtmIndx[i] = qmAtmIndxVec[i] ;
510 std::map<int, std::vector<std::pair<int,int> > > qmGrpIDBonds ;
518 for (
int bndIt = 0; bndIt <
numBonds; bndIt++) {
520 bond curr = bonds[bndIt] ;
523 if ( qmBondValVec[curr.
atom1] != 0 &&
524 qmBondValVec[curr.
atom2] != 0 ) {
527 if (qmAtomGroup[curr.
atom1] != 0 &&
528 qmAtomGroup[curr.
atom2] != 0) {
530 curr.
atom2 <<
" are assigned as QM atoms.\n" <<
endi;
531 NAMD_die(
"Error in QM-MM bond assignment.") ;
535 if (qmAtomGroup[curr.
atom1] == 0 &&
536 qmAtomGroup[curr.
atom2] == 0) {
538 curr.
atom2 <<
" are assigned as MM atoms.\n" <<
endi;
539 NAMD_die(
"Error in QM-MM bond assignment.") ;
543 std::pair<int,int> newPair(0,0);
546 if (qmAtomGroup[curr.
atom1] != 0) {
547 newPair = std::pair<int,int>(curr.
atom2, curr.
atom1) ;
548 currGrpID = qmAtomGroup[curr.
atom1];
550 newPair = std::pair<int,int>(curr.
atom1, curr.
atom2) ;
551 currGrpID = qmAtomGroup[curr.
atom2];
554 int grpIndx = qmGrpIDMap[currGrpID] ;
557 auto retIter = qmGrpIDBonds.find(grpIndx) ;
560 if (retIter == qmGrpIDBonds.end()) {
561 qmGrpIDBonds.insert(std::pair<
BigReal,std::vector<std::pair<int,int> > >(
562 grpIndx, std::vector<std::pair<int,int> >() ) );
565 qmGrpIDBonds[grpIndx].push_back( newPair );
575 if (bondCounter == 0)
576 iout <<
iWARN <<
"We found " << bondCounter <<
" QM-MM bonds.\n" <<
endi ;
578 iout <<
iINFO <<
"We found " << bondCounter <<
" QM-MM bonds.\n" <<
endi ;
583 qmNumBonds = bondCounter ;
585 qmMMBond =
new int*[qmNumBonds];
587 qmDummyBondVal =
new BigReal[qmNumBonds];
589 qmMMChargeTarget =
new int*[qmNumBonds] ;
590 qmMMNumTargs =
new int[qmNumBonds] ;
592 qmDummyElement =
new char*[qmNumBonds] ;
595 qmNumGrps = atmGrpSet.size();
597 qmGrpSizes =
new int[qmNumGrps];
599 qmGrpID =
new Real[qmNumGrps];
601 qmGrpChrg =
new Real[qmNumGrps];
603 qmGrpMult =
new Real[qmNumGrps] ;
605 qmGrpNumBonds =
new int[qmNumGrps];
607 qmGrpBonds =
new int*[qmNumGrps];
608 qmMMBondedIndx =
new int*[qmNumGrps];
617 for (
int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
618 qmGrpMult[grpIter] = 0;
623 for ( ; current; current = current->
next ) {
625 auto strVec =
split( std::string(current->
data) ,
" ");
627 if (strVec.size() != 2 ) {
628 iout <<
iERROR <<
"Format error in QM Multiplicity string: "
631 NAMD_die(
"Error processing QM information.");
634 std::stringstream storConv ;
636 storConv << strVec[0] ;
639 if (storConv.fail()) {
640 iout <<
iERROR <<
"Error parsing QMMult selection: "
643 NAMD_die(
"Error processing QM information.");
647 storConv << strVec[1];
649 storConv >> multiplicity;
650 if (storConv.fail()) {
651 iout <<
iERROR <<
"Error parsing QMMult selection: "
654 NAMD_die(
"Error processing QM information.");
657 auto it = qmGrpIDMap.find(grpID);
659 if (it == qmGrpIDMap.end()) {
660 iout <<
iERROR <<
"Error parsing QMMult selection. Could not find QM group ID: "
663 NAMD_die(
"Error processing QM information.");
666 iout <<
iINFO <<
"Applying user defined multiplicity "
667 << multiplicity <<
" to QM group ID " << grpID <<
"\n" <<
endi;
668 qmGrpMult[it->second] = multiplicity;
675 NAMD_die(
"ORCA neds multiplicity values for all QM regions.");
683 for (
size_t grpIndx = 0; grpIndx < qmGrpSizeVec.size(); grpIndx++) {
685 bool nonInteger =
true;
686 if ((fabsf(roundf(grpChrgVec[grpIndx]) - grpChrgVec[grpIndx]) <= 0.001f) ) {
687 grpChrgVec[grpIndx] = roundf(grpChrgVec[grpIndx]) ;
691 iout <<
iINFO << grpIndx + 1 <<
") Group ID: " << qmGroupIDsVec[grpIndx]
692 <<
" ; Group size: " << qmGrpSizeVec[grpIndx] <<
" atoms"
693 <<
" ; Total PSF charge: " << grpChrgVec[grpIndx] <<
"\n" <<
endi ;
695 if (nonInteger && simParams->
PMEOn)
696 NAMD_die(
"QM atoms do not add up to a whole charge, which is needed for PME.") ;
700 current = cfgList->
find(
"QMCharge");
701 for ( ; current; current = current->
next ) {
703 auto strVec =
split( std::string(current->
data) ,
" ");
705 if (strVec.size() != 2 ) {
706 iout <<
iERROR <<
"Format error in QM Charge string: "
709 NAMD_die(
"Error processing QM information.");
712 std::stringstream storConv ;
714 storConv << strVec[0] ;
717 if (storConv.fail()) {
718 iout <<
iERROR <<
"Error parsing QMCharge selection: "
721 NAMD_die(
"Error processing QM information.");
725 storConv << strVec[1];
728 if (storConv.fail()) {
729 iout <<
iERROR <<
"Error parsing QMCharge selection: "
732 NAMD_die(
"Error processing QM information.");
735 auto it = qmGrpIDMap.find(grpID);
737 if (it == qmGrpIDMap.end()) {
738 iout <<
iERROR <<
"Error parsing QMCharge selection. Could not find QM group ID: "
741 NAMD_die(
"Error processing QM information.");
744 iout <<
iINFO <<
"Found user defined charge "
745 << charge <<
" for QM group ID " << grpID <<
". Will ignore PSF charge.\n" <<
endi;
746 grpChrgVec[it->second] =
charge;
758 std::string::size_type chrgLoc = std::string::npos ;
761 std::string configLine(cfgList->
find(
"QMConfigLine")->
data) ;
762 chrgLoc = configLine.find(
"CHARGE") ;
764 if ( chrgLoc != std::string::npos ) {
765 iout <<
iINFO <<
"Found user defined charge in command line. This \
766 will be used for all QM regions and will take precedence over all other charge \
767 definitions.\n" <<
endi;
769 else if ( chrgLoc == std::string::npos && (simParams->
qmChrgFromPSF || chrgCount == qmNumGrps) ) {
780 iout <<
iERROR <<
"Could not find charge for all QM groups. For MOPAC, \
781 charges can be defined through QMConfigLine, QMCharge or QMChargeFromPSF keywords.\n" <<
endi;
782 NAMD_die(
"Error processing QM information.");
788 if ((chrgCount != qmNumGrps ) && !simParams->
qmChrgFromPSF) {
792 iout <<
iERROR <<
"Could not find charge for all QM groups. For ORCA, \
793 charges can be defined through QMCharge or QMChargeFromPSF keywords.\n" <<
endi;
794 NAMD_die(
"Error processing QM information.");
800 if (qmNumBonds > 0 && qmNoPC) {
801 qmMeNumBonds = qmNumBonds;
802 qmMeMMindx =
new int[qmMeNumBonds] ;
803 qmMeQMGrp =
new Real[qmMeNumBonds] ;
817 for (
int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
819 qmGrpID[grpIter] = qmGroupIDsVec[grpIter] ;
821 qmGrpSizes[grpIter] = qmGrpSizeVec[grpIter] ;
823 qmGrpChrg[grpIter] = grpChrgVec[grpIter];
827 int currNumbBonds = qmGrpIDBonds[grpIter].size() ;
830 qmGrpNumBonds[grpIter] = currNumbBonds;
832 if (currNumbBonds > 0) {
834 qmGrpBonds[grpIter] =
new int[currNumbBonds];
835 qmMMBondedIndx[grpIter] =
new int[currNumbBonds];
837 for (
int bndIter = 0; bndIter<currNumbBonds; bndIter++) {
840 qmMMBond[bondCounter] =
new int[2] ;
841 qmMMBond[bondCounter][0] = qmGrpIDBonds[grpIter][bndIter].first ;
842 qmMMBond[bondCounter][1] = qmGrpIDBonds[grpIter][bndIter].second ;
845 qmGrpBonds[grpIter][bndIter] = bondCounter;
848 qmMMBondedIndx[grpIter][bndIter] = qmGrpIDBonds[grpIter][bndIter].first ;
851 qmDummyElement[bondCounter] =
new char[3];
852 strcpy(qmDummyElement[bondCounter],
"H\0");
858 if (qmBondValVec[qmMMBond[bondCounter][0]] !=
859 qmBondValVec[qmMMBond[bondCounter][1]]
861 iout <<
iERROR <<
"qmBondDist is ON but the values in the bond column are different for atom "
862 << qmMMBond[bondCounter][0] <<
" and " << qmMMBond[bondCounter][1]
864 NAMD_die(
"Error in QM-MM bond processing.");
867 bondVal = qmBondValVec[qmMMBond[bondCounter][0]] ;
871 if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],
"C") == 0 ) {
874 else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],
"O") == 0 ) {
877 else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],
"N") == 0 ) {
881 iout <<
iERROR <<
"qmBondDist is OFF but the QM atom has no default distance value. Atom "
882 << qmMMBond[bondCounter][1] <<
", with element: "
883 << qmElementArray[qmMMBond[bondCounter][1]]
885 NAMD_die(
"Error in QM-MM bond processing.");
890 iout <<
iINFO <<
"MM-QM pair: " << qmMMBond[bondCounter][0] <<
":"
891 << qmMMBond[bondCounter][1]
892 <<
" -> Value (distance or ratio): " << bondVal
893 <<
" (QM Group " << grpIter <<
" ID " << qmGrpID[grpIter] <<
")"
896 qmDummyBondVal[bondCounter] = bondVal;
900 if (qmMeNumBonds > 0) {
901 qmMeMMindx[bondCounter] = qmMMBond[bondCounter][0];
902 qmMeQMGrp[bondCounter] = qmGrpID[grpIter];
917 current = cfgList->
find(
"QMLinkElement");
919 int numParsedLinkElem = 0;
920 for ( ; current != NULL; current = current->
next ) {
922 DebugM(3,
"Parsing link atom element: " << current->
data <<
"\n" );
924 strVec =
split( std::string(current->
data) ,
" ");
928 if (strVec.size() != 3 && qmNumBonds > 1) {
929 iout <<
iERROR <<
"Format error in QM link atom element selection: "
932 NAMD_die(
"Error processing QM information.");
936 if (strVec.size() != 1 && strVec.size() != 3 && qmNumBonds == 1) {
937 iout <<
iERROR <<
"Format error in QM link atom element selection: "
940 NAMD_die(
"Error processing QM information.");
943 std::stringstream storConv ;
944 int bondAtom1, bondAtom2;
945 std::string linkElement ;
947 if (strVec.size() == 1) {
948 linkElement = strVec[0].substr(0,2);
950 else if (strVec.size() == 3) {
952 storConv << strVec[0] ;
953 storConv >> bondAtom1;
954 if (storConv.fail()) {
955 iout <<
iERROR <<
"Error parsing link atom element selection: "
958 NAMD_die(
"Error processing QM information.");
962 storConv << strVec[1];
963 storConv >> bondAtom2;
964 if (storConv.fail()) {
965 iout <<
iERROR <<
"Error parsing link atom element selection: "
968 NAMD_die(
"Error processing QM information.");
971 linkElement = strVec[2].substr(0,2);
976 if (numParsedLinkElem>qmNumBonds) {
977 iout <<
iERROR <<
"More link atom elements were set than there"
978 " are QM-MM bonds.\n" <<
endi;
979 NAMD_die(
"Error processing QM information.");
984 if (strVec.size() == 1) {
987 else if (strVec.size() == 3) {
989 Bool foundBond =
false;
991 for (bondIter=0; bondIter<qmNumBonds; bondIter++) {
993 if ( (qmMMBond[bondIter][0] == bondAtom1 &&
994 qmMMBond[bondIter][1] == bondAtom2 ) ||
995 (qmMMBond[bondIter][0] == bondAtom2 &&
996 qmMMBond[bondIter][1] == bondAtom1 ) ) {
1004 iout <<
iERROR <<
"Error parsing link atom element selection: "
1007 iout <<
iERROR <<
"Atom pair was not found to be a QM-MM bond.\n"
1009 NAMD_die(
"Error processing QM information.");
1013 strcpy(qmDummyElement[bondIter],linkElement.c_str());
1014 qmDummyElement[bondIter][2] =
'\0';
1016 iout <<
iINFO <<
"Applying user defined link atom element "
1017 << qmDummyElement[bondIter] <<
" to QM-MM bond "
1018 << bondIter <<
": " << qmMMBond[bondIter][0]
1019 <<
" - " << qmMMBond[bondIter][1]
1029 int32 **bondsWithAtomLocal = NULL ;
1030 int32 *byAtomSizeLocal = NULL;
1032 if (qmNumBonds > 0) {
1041 byAtomSizeLocal[i] = 0;
1045 byAtomSizeLocal[bonds[i].
atom1]++;
1046 byAtomSizeLocal[bonds[i].atom2]++;
1050 bondsWithAtomLocal[i] = tmpArenaLocal->
getNewArray(byAtomSizeLocal[i]+1);
1051 bondsWithAtomLocal[i][byAtomSizeLocal[i]] = -1;
1052 byAtomSizeLocal[i] = 0;
1056 int a1 = bonds[i].
atom1;
1057 int a2 = bonds[i].
atom2;
1058 bondsWithAtomLocal[a1][byAtomSizeLocal[a1]++] = i;
1059 bondsWithAtomLocal[a2][byAtomSizeLocal[a2]++] = i;
1066 for (
int qmBndIt = 0; qmBndIt < qmNumBonds; qmBndIt++) {
1070 std::vector<int> chrgTrgt ;
1072 int MM1 = qmMMBond[qmBndIt][0], MM2, MM2BondIndx, MM3, MM3BondIndx;
1081 for (
int i=0; i<byAtomSizeLocal[MM1]; i++) {
1082 MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1086 if (bonds[MM2BondIndx].atom1 == MM1)
1087 MM2 = bonds[MM2BondIndx].
atom2;
1089 MM2 = bonds[MM2BondIndx].
atom1;
1093 if (qmAtomGroup[MM2] > 0)
1096 chrgTrgt.push_back(MM2);
1104 for (
int i=0; i<byAtomSizeLocal[MM1]; i++) {
1105 MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1109 if (bonds[MM2BondIndx].atom1 == MM1)
1110 MM2 = bonds[MM2BondIndx].
atom2;
1112 MM2 = bonds[MM2BondIndx].
atom1;
1116 if (qmAtomGroup[MM2] > 0)
1119 for (
int i=0; i<byAtomSizeLocal[MM2]; i++) {
1120 MM3BondIndx = bondsWithAtomLocal[MM2][i] ;
1124 if (bonds[MM3BondIndx].atom1 == MM2)
1125 MM3 = bonds[MM3BondIndx].
atom2;
1127 MM3 = bonds[MM3BondIndx].
atom1;
1132 if (qmAtomGroup[MM3] > 0 || MM3 == MM1)
1135 chrgTrgt.push_back(MM3);
1145 for (
int i=0; i<byAtomSizeLocal[MM1]; i++) {
1146 MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1150 if (bonds[MM2BondIndx].atom1 == MM1)
1151 MM2 = bonds[MM2BondIndx].
atom2;
1153 MM2 = bonds[MM2BondIndx].
atom1;
1157 if (qmAtomGroup[MM2] > 0)
1160 chrgTrgt.push_back(MM2);
1168 chrgTrgt.push_back(MM1);
1173 qmMMChargeTarget[qmBndIt] =
new int[chrgTrgt.size()] ;
1174 qmMMNumTargs[qmBndIt] = chrgTrgt.size();
1176 DebugM(3,
"MM-QM bond " << qmBndIt <<
"; MM atom "
1177 << qmMMBond[qmBndIt][0] <<
" conections: \n" );
1179 for (
size_t i=0; i < chrgTrgt.size(); i++) {
1180 qmMMChargeTarget[qmBndIt][i] = chrgTrgt[i];
1181 DebugM(3,
"MM Bonded to: " << chrgTrgt[i] <<
"\n" );
1187 if (bondsWithAtomLocal != NULL)
1188 delete [] bondsWithAtomLocal; bondsWithAtomLocal = 0;
1189 if (byAtomSizeLocal != NULL)
1190 delete [] byAtomSizeLocal; byAtomSizeLocal = 0;
1191 if (tmpArenaLocal != NULL)
1192 delete tmpArenaLocal; tmpArenaLocal = 0;
1201 std::map<Real,int> grpLSSSize ;
1202 std::map<Real,int>::iterator itGrpSize;
1204 qmLSSTotalNumAtms = 0;
1205 qmLSSResidueSize = 0;
1216 std::map<Real, int> grpNumLSSRes;
1217 std::map<Real, int>::iterator itGrpNumRes;
1219 for(
auto it = lssGrpRes.
begin(); it != lssGrpRes.
end();it++ ) {
1221 if (it->atmIDs.size() != it->size) {
1222 iout <<
iERROR <<
"The number of atoms loaded for residue "
1223 << it->resID <<
" does not match the expected for this residue type.\n"
1225 NAMD_die(
"Error parsing data for LSS.");
1228 qmLSSTotalNumAtms += it->size;
1231 if (resSize < 0) resSize = it->size ;
1232 if (resSize > 0 and resSize != it->size) {
1233 iout <<
iERROR <<
"The number of atoms loaded for residue "
1234 << it->resID <<
" does not match previously loaded residues.\n"
1236 NAMD_die(
"Error parsing data for LSS.");
1248 itGrpSize = grpLSSSize.find(it->qmGrpID) ;
1249 if (itGrpSize != grpLSSSize.end())
1250 itGrpSize->second += it->size;
1252 grpLSSSize.insert(std::pair<Real,int>(it->qmGrpID, it->size));
1255 itGrpNumRes = grpNumLSSRes.find(it->qmGrpID) ;
1256 if (itGrpNumRes != grpNumLSSRes.end())
1257 itGrpNumRes->second += 1;
1259 grpNumLSSRes.insert(std::pair<Real,int>(it->qmGrpID, 1));
1262 qmLSSResidueSize = lssGrpRes.
begin()->size ;
1264 qmLSSSize =
new int[qmNumGrps];
1266 qmLSSIdxs =
new int[qmLSSTotalNumAtms];
1274 qmLSSRefSize =
new int[qmNumGrps];
1276 qmLSSMass =
new Mass[qmLSSTotalNumAtms];
1278 qmLSSRefIDs =
new int[totalNumRefAtms];
1279 qmLSSRefMass =
new Mass[totalNumRefAtms];
1282 for (
int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
1284 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
1286 if (itGrpSize != grpNumLSSRes.end())
1287 qmLSSSize[grpIndx] = itGrpSize->second;
1289 qmLSSSize[grpIndx] = 0;
1291 for(
auto it = lssGrpRes.
begin(); it != lssGrpRes.
end();it++ ) {
1293 if (it->qmGrpID == qmGrpID[grpIndx]) {
1294 for (
int i=0; i<it->atmIDs.size(); i++) {
1295 qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
1296 qmLSSMass[lssAtmIndx] = atoms[it->atmIDs[i]].
mass;
1302 DebugM(4,
"QM group " << qmGrpID[grpIndx]
1303 <<
" has " << qmLSSSize[grpIndx] <<
" solvent molecules, "
1304 <<
"totalling " << grpLSSSize[qmGrpID[grpIndx]]
1305 <<
" atoms (check " << lssAtmIndx <<
").\n" );
1307 qmLSSRefSize[grpIndx] = lssGrpRefIDs[qmGrpID[grpIndx]].size();
1308 for(
int i=0; i < qmLSSRefSize[grpIndx]; i++) {
1309 qmLSSRefIDs[lssRefIndx] = lssGrpRefIDs[qmGrpID[grpIndx]][i];
1310 qmLSSRefMass[lssRefIndx] = atoms[qmLSSRefIDs[lssRefIndx]].
mass;
1314 DebugM(4,
"QM group " << qmGrpID[grpIndx]
1315 <<
" has " << qmLSSRefSize[grpIndx] <<
" non-solvent atoms for LSS.\n" );
1322 for (
int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
1324 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
1326 if (itGrpSize != grpNumLSSRes.end())
1327 qmLSSSize[grpIndx] = itGrpSize->second;
1329 qmLSSSize[grpIndx] = 0;
1331 for(
auto it = lssGrpRes.
begin(); it != lssGrpRes.
end();it++ ) {
1333 if (it->qmGrpID == qmGrpID[grpIndx]) {
1334 for (
int i=0; i<it->atmIDs.size(); i++) {
1335 qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
1341 DebugM(4,
"QM group " << qmGrpID[grpIndx]
1342 <<
" has " << qmLSSSize[grpIndx] <<
" solvent molecules, "
1343 <<
"totalling " << grpLSSSize[qmGrpID[grpIndx]]
1344 <<
" atoms (check " << lssAtmIndx <<
").\n" );
1363 current = cfgList->
find(
"QMCustomPCFile");
1366 std::map<Real,std::vector<int> > qmPCVecMap ;
1368 int numParsedPBDs = 0;
1369 for ( ; current != NULL; current = current->
next ) {
1372 customPCPDB =
new PDB(current->
data);
1375 NAMD_die(
"Number of atoms in QM Custom PC PDB file doesn't match coordinate PDB");
1377 std::vector< int > currPCSel ;
1379 int currGrpSize = 0 ;
1381 for (
int atmInd = 0 ; atmInd <
numAtoms; atmInd++) {
1386 if ( beta != 0 && occ != 0)
1387 NAMD_die(
"An atom cannot be marked as QM and poitn charge simultaneously!");
1391 currPCSel.push_back(atmInd) ;
1396 NAMD_die(
"QM Group selection is different in reference PDB and Custom-PC PDB!");
1398 if (currQMID == 0) {
1404 if (currQMID != beta)
1405 NAMD_die(
"Found two different QM group IDs in this file!");
1413 if (currGrpSize != qmGrpSizeVec[ qmGrpIDMap[currQMID] ])
1414 NAMD_die(
"Reference PDB and Custom-PC PDB have different QM group sizes!") ;
1416 qmPCVecMap.insert(std::pair<
Real,std::vector<int> >(
1417 currQMID, currPCSel ));
1425 if (numParsedPBDs != qmNumGrps && simParams->
qmCustomPCSel) {
1426 iout <<
iWARN <<
"The number of files provided for custom point "
1427 "charges is not equal to the number of QM groups!\n" <<
endi;
1432 qmCustPCSizes =
new int[qmNumGrps];
1433 for (
int i=0; i<qmNumGrps; i++)
1434 qmCustPCSizes[i] = 0;
1440 for (
auto it = qmPCVecMap.begin(); it != qmPCVecMap.end(); it++) {
1441 qmTotCustPCs += (*it).second.size();
1442 int qmIndx = qmGrpIDMap[(*it).first];
1443 qmCustPCSizes[qmIndx] = (*it).second.size();
1446 qmCustomPCIdxs =
new int[qmTotCustPCs];
1451 for (
int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
1453 Real currQMID = qmGrpID[grpIndx];
1455 iout <<
iINFO <<
"Loading " << qmPCVecMap[currQMID].size()
1456 <<
" custom point charges to QM Group " << grpIndx
1457 <<
" (ID: " << currQMID <<
")\n" <<
endi;
1459 for (
int pcIndx=0; pcIndx<qmPCVecMap[currQMID].size(); pcIndx++) {
1460 qmCustomPCIdxs[auxIter] = qmPCVecMap[currQMID][pcIndx] ;
1467 qmcSMD = simParams->
qmCSMD;
1469 read_qm_csdm_file(qmGrpIDMap);
1481 for (
int i=0; i<qmNumQMAtoms; i++) {
1482 qmAtmChrg[i] = atoms[qmAtmIndx[i]].
charge;
1483 atoms[qmAtmIndx[i]].
charge = 0;
1494 int a1,a2;
float k, ref;
1496 for ( ; file; file = file->
next ) {
1497 FILE *f = fopen(file->
data,
"r");
1511 if (ret_code!=0)
break;
1514 sscanf(buffer,
"%s",type);
1517 if ( ! strncasecmp(type,
"bond",4) ) {
1518 if ( sscanf(buffer,
"%s %d %d %f %f %s",
1519 type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1;
1524 if( qmAtomGroup[a1] > 0 && qmAtomGroup[a2]) {
1528 qmExtraBonds.add(tmp);
1532 }
else if ( ! strncasecmp(type,
"#",1) ) {
1543 #endif // MEM_OPT_VERSION
1549 #ifdef MEM_OPT_VERSION
1550 NAMD_die(
"QMMM interface is not supported in memory optimized builds.");
1553 DebugM(3,
"Cleaning QM bonds, angles, dihedrals, impropers and crossterms.\n")
1558 int nonQMBondCount = 0;
1560 for (
int i = 0; i <
numBonds; i++) {
1562 int part1 = qmAtomGroup[bonds[i].
atom1];
1563 int part2 = qmAtomGroup[bonds[i].atom2];
1564 if (part1 > 0 && part2 > 0 ) {
1571 for (
int ebi=0; ebi < qmExtraBonds.size() ; ebi++) {
1573 if( (qmExtraBonds[ebi].atom1 == bonds[i].atom1 ||
1574 qmExtraBonds[ebi].atom1 == bonds[i].atom2) &&
1575 (qmExtraBonds[ebi].atom2 == bonds[i].atom1 ||
1576 qmExtraBonds[ebi].atom2 == bonds[i].atom2) ) {
1578 nonQMBonds[nonQMBondCount++] = bonds[i];
1589 if ((part1 > 0 || part2 > 0) && qmNumBonds == 0) {
1590 iout <<
iERROR <<
" Atoms " << bonds[i].atom1
1591 <<
" and " << bonds[i].atom2 <<
" form a QM-MM bond!\n" <<
endi ;
1592 NAMD_die(
"No QM-MM bonds were defined by the user, but one was found.");
1594 nonQMBonds[nonQMBondCount++] = bonds[i];
1597 numBonds = nonQMBondCount;
1600 for (
int i = 0; i < nonQMBondCount; i++) {
1601 bonds[i]=nonQMBonds[i];
1603 delete [] nonQMBonds;
1609 int nonQMAngleCount = 0;
1610 qmDroppedAngles = 0;
1612 int part1 = qmAtomGroup[angles[i].
atom1];
1613 int part2 = qmAtomGroup[angles[i].atom2];
1614 int part3 = qmAtomGroup[angles[i].atom3];
1615 if (part1 > 0 && part2 > 0 && part3 > 0) {
1620 nonQMAngles[nonQMAngleCount++] = angles[i];
1623 numAngles = nonQMAngleCount;
1626 for (
int i = 0; i < nonQMAngleCount; i++) {
1627 angles[i]=nonQMAngles[i];
1629 delete [] nonQMAngles;
1635 int nonQMDihedralCount = 0;
1636 qmDroppedDihedrals = 0;
1638 int part1 = qmAtomGroup[dihedrals[i].
atom1];
1639 int part2 = qmAtomGroup[dihedrals[i].atom2];
1640 int part3 = qmAtomGroup[dihedrals[i].atom3];
1641 int part4 = qmAtomGroup[dihedrals[i].atom4];
1642 if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
1644 qmDroppedDihedrals++;
1647 nonQMDihedrals[nonQMDihedralCount++] = dihedrals[i];
1650 numDihedrals = nonQMDihedralCount;
1651 delete [] dihedrals;
1654 dihedrals[i]=nonQMDihedrals[i];
1656 delete [] nonQMDihedrals;
1661 int nonQMImproperCount = 0;
1662 qmDroppedImpropers = 0;
1664 int part1 = qmAtomGroup[impropers[i].
atom1];
1665 int part2 = qmAtomGroup[impropers[i].atom2];
1666 int part3 = qmAtomGroup[impropers[i].atom3];
1667 int part4 = qmAtomGroup[impropers[i].atom4];
1668 if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
1670 qmDroppedImpropers++;
1673 nonQMImpropers[nonQMImproperCount++] = impropers[i];
1676 numImpropers = nonQMImproperCount;
1677 delete [] impropers;
1680 impropers[i]=nonQMImpropers[i];
1682 delete [] nonQMImpropers;
1687 int nonQMCrosstermCount = 0;
1688 qmDroppedCrossterms = 0 ;
1690 int part1 = qmAtomGroup[crossterms[i].
atom1];
1691 int part2 = qmAtomGroup[crossterms[i].atom2];
1692 int part3 = qmAtomGroup[crossterms[i].atom3];
1693 int part4 = qmAtomGroup[crossterms[i].atom4];
1694 int part5 = qmAtomGroup[crossterms[i].atom5];
1695 int part6 = qmAtomGroup[crossterms[i].atom6];
1696 int part7 = qmAtomGroup[crossterms[i].atom7];
1697 int part8 = qmAtomGroup[crossterms[i].atom8];
1698 if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0 &&
1699 part5 > 0 && part6 > 0 && part7 > 0 && part8 > 0) {
1701 qmDroppedCrossterms++;
1704 nonQMCrossterms[nonQMCrosstermCount++] = crossterms[i];
1707 numCrossterms = nonQMCrosstermCount;
1708 delete [] crossterms;
1711 crossterms[i] = nonQMCrossterms[i] ;
1713 delete [] nonQMCrossterms ;
1716 iout <<
iINFO <<
"The QM region will remove " << qmDroppedBonds <<
" bonds, " <<
1717 qmDroppedAngles <<
" angles, " << qmDroppedDihedrals <<
" dihedrals, "
1718 << qmDroppedImpropers <<
" impropers and " << qmDroppedCrossterms
1719 <<
" crossterms.\n" <<
endi ;
1727 void Molecule::read_qm_csdm_file(std::map<Real,int> &qmGrpIDMap) {
1729 std::ifstream cSMDInputFile ;
1730 std::string Line(
"") ;
1732 iout <<
iINFO <<
"Reading QM cSMD configuration file: " <<
1735 cSMDInputFile.open( simParams->
qmCSMDFile ) ;
1737 if (! cSMDInputFile.good())
1739 NAMD_die(
"Configuration file for QM cSMD could not be read!") ;
1745 cSMDindexV.
resize(qmNumGrps);
1755 while( ! cSMDInputFile.eof() )
1757 getline(cSMDInputFile, Line) ;
1759 if ( ! Line.length() )
1762 if (Line.substr(0,1) == std::string(
"#") )
1765 auto strVec =
split( Line,
" ");
1767 if (strVec.size() != 5 ) {
1768 iout <<
iERROR <<
"Format error in QM cSDM configuration file: "
1769 << Line <<
"\n" <<
endi;
1770 NAMD_die(
"Error processing QM information.");
1773 std::stringstream storConv ;
1781 for (
int i=0; i < 5; i++ ) {
1783 storConv << strVec[i] ;
1784 storConv >> convData;
1785 if (storConv.fail()) {
1786 iout <<
iERROR <<
"Error parsing QM cSMD configuration file: "
1787 << convData <<
"\n" <<
endi;
1788 NAMD_die(
"Error processing QM information.");
1793 cSMDpairsV[cSMDnumInst].first = convData ;
1796 cSMDpairsV[cSMDnumInst].second = convData ;
1799 cSMDKsV[cSMDnumInst] = convData ;
1802 cSMDVelsV[cSMDnumInst] = convData ;
1805 cSMDcoffsV[cSMDnumInst] = convData ;
1811 if (cSMDpairsV[cSMDnumInst].first == cSMDpairsV[cSMDnumInst].second) {
1812 iout <<
iERROR <<
"Conditional SMD atoms must be different! We got " <<
1813 cSMDpairsV[cSMDnumInst].first <<
" and " << cSMDpairsV[cSMDnumInst].second <<
"\n" <<
endi;
1814 NAMD_die(
"Error processing QM information.") ;
1818 if (qmAtomGroup[cSMDpairsV[cSMDnumInst].first] == 0 ||
1819 qmAtomGroup[cSMDpairsV[cSMDnumInst].second] == 0) {
1820 iout <<
iERROR <<
"Atoms " << cSMDpairsV[cSMDnumInst].first <<
" and " <<
1821 cSMDpairsV[cSMDnumInst].second <<
" MUST be assigned as QM atoms.\n" <<
endi;
1822 NAMD_die(
"Error processing QM information.") ;
1826 if (qmAtomGroup[cSMDpairsV[cSMDnumInst].first] !=
1827 qmAtomGroup[cSMDpairsV[cSMDnumInst].second] ) {
1828 iout <<
iERROR <<
"Atoms in cSMD MUST be assigned to the same QM group.\n" <<
endi;
1829 NAMD_die(
"Error processing QM information.") ;
1832 int grpIndx = qmGrpIDMap[qmAtomGroup[cSMDpairsV[cSMDnumInst].first]] ;
1834 iout <<
iINFO <<
"Adding cSMD data: (" << cSMDnumInst <<
") " <<
1835 cSMDpairsV[cSMDnumInst].first <<
"," << cSMDpairsV[cSMDnumInst].second <<
"," <<
1836 cSMDKsV[cSMDnumInst] <<
"," << cSMDVelsV[cSMDnumInst] <<
"," << cSMDcoffsV[cSMDnumInst] <<
1837 " to QM Group " << grpIndx <<
"\n" <<
endi ;
1839 cSMDindexV[grpIndx].
add(cSMDnumInst);
1844 cSMDindex =
new int*[qmNumGrps];
1845 cSMDindxLen =
new int[qmNumGrps];
1847 for (
int i=0; i<qmNumGrps; i++) {
1848 cSMDindex[i] =
new int[cSMDindexV[i].
size()];
1849 cSMDindxLen[i] = cSMDindexV[i].size();
1851 for (
int j=0; j<cSMDindxLen[i]; j++)
1852 cSMDindex[i][j] = cSMDindexV[i][j] ;
1855 cSMDpairs =
new int*[cSMDnumInst];
1856 for (
int i=0; i<cSMDnumInst; i++)
1857 cSMDpairs[i] =
new int[2];
1858 cSMDKs =
new Real[cSMDnumInst];
1859 cSMDVels =
new Real[cSMDnumInst];
1860 cSMDcoffs =
new Real[cSMDnumInst];
1862 for (
int i=0; i<cSMDnumInst; i++) {
1863 cSMDpairs[i][0] = cSMDpairsV[i].first;
1864 cSMDpairs[i][1] = cSMDpairsV[i].second;
1866 cSMDKs[i] = cSMDKsV[i];
1867 cSMDVels[i] = cSMDVelsV[i];
1868 cSMDcoffs[i] = cSMDcoffsV[i];
int get_residue_size(const char *segid, int resid) const
std::ostream & iINFO(std::ostream &s)
Type * getNewArray(int n)
std::map< Real, refSelStrVec > refSelStrMap
int NAMD_read_line(FILE *fd, char *buf, int bufsize)
void assign_vdw_index(const char *, Atom *)
std::vector< int > atmIDs
BigReal temperaturefactor(void)
std::ostream & endi(std::ostream &s)
char qmCSMDFile[NAMD_FILENAME_BUFFER_SIZE]
std::ostream & iWARN(std::ostream &s)
refSelStr(std::string newSegID, int newResid)
int insert(const Elem &elem)
const char * residuename(void)
PDBAtom * atom(int place)
int NAMD_blank_string(char *str)
std::pair< Real, refSelStrVec > refSelStrPair
int lookup(const char *segid, int resid, int *begin, int *end) const
void prepare_qm(const char *pdbFileName, Parameters *params, ConfigList *cfgList)
void NAMD_die(const char *err_msg)
bool operator<(const qmSolvData &ref)
std::vector< std::string > split(const std::string &text, std::string delimiter)
int add(const Elem &elem)
Bool qmMOPACAddConfigChrg
std::vector< refSelStr > refSelStrVec
bool operator==(const qmSolvData &ref)
StringList * find(const char *name) const
Elem * find(const Elem &elem)
std::pair< Origin, Target > cSMDPair
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
std::ostream & iERROR(std::ostream &s)
qmSolvData(int newResID, int newBegAtm, int newSize, const char *newSegName, Real newQMID)
const char * element(void)
const char * segmentname(void)
qmSolvData(int newResID, int newBegAtm, int newSize)