MoleculeQM.C

Go to the documentation of this file.
00001 
00007 #ifdef DEBUG_QM
00008   #define DEBUGM
00009 #endif
00010 
00011 #include <stdio.h>
00012 #include <string.h>
00013 #include <stdlib.h>
00014 #include <ctype.h>
00015 #include "ResizeArray.h"
00016 #include "InfoStream.h"
00017 #include "Molecule.h"
00018 #include "strlib.h"
00019 #include "MStream.h"
00020 #include "Communicate.h"
00021 #include "Node.h"
00022 #include "ObjectArena.h"
00023 #include "Parameters.h"
00024 #include "PDB.h"
00025 #include "SimParameters.h"
00026 #include "Hydrogen.h"
00027 #include "UniqueSetIter.h"
00028 #include "charm++.h"
00029 
00030 #define MIN_DEBUG_LEVEL 3
00031 //#define DEBUGM
00032 #include "Debug.h"
00033 #include "CompressPsf.h"
00034 #include "ParallelIOMgr.h"
00035 #include <deque>
00036 #include <algorithm>
00037 
00038 #ifndef CODE_REDUNDANT
00039 #define CODE_REDUNDANT 0
00040 #endif
00041 
00042 #include <string>
00043 #include <sstream>
00044 #include <fstream>
00045 
00046 class qmSolvData {
00047 public:
00048     char segName[5];
00049     int resID, begAtmID, size;
00050     std::vector<int> atmIDs ;
00051     Real qmGrpID ;
00052     
00053     qmSolvData() : resID(-1), begAtmID(-1), size(0) {}
00054     qmSolvData(int newResID, int newBegAtm, int newSize) {
00055         resID = newResID;
00056         begAtmID = newBegAtm;
00057         size = newSize;
00058     }
00059     qmSolvData(int newResID, int newBegAtm, int newSize, 
00060                const char* newSegName, Real newQMID) {
00061         resID = newResID;
00062         begAtmID = newBegAtm;
00063         size = newSize;
00064         strncpy(segName, newSegName,5);
00065         atmIDs.push_back(newBegAtm) ;
00066         qmGrpID = newQMID;
00067     }
00068     
00069     bool operator <(const qmSolvData& ref) {return begAtmID < ref.begAtmID;}
00070     bool operator ==(const qmSolvData& ref) {return begAtmID == ref.begAtmID;}
00071 };
00072 
00073 std::vector<std::string> split(const std::string &text, std::string delimiter) {
00074     
00075     std::vector<std::string> tokens;
00076     std::size_t start = 0, end = 0;
00077     
00078     while ((end = text.find(delimiter, start)) != std::string::npos) {
00079         
00080         std::string temp = text.substr(start, end - start);
00081         
00082         if (! temp.empty()) tokens.push_back(temp);
00083         
00084         start = end + delimiter.length();
00085     }
00086     
00087     // Gets last item
00088     std::string temp = text.substr(start);
00089     if (! temp.empty()) tokens.push_back(temp);
00090     
00091     return tokens;
00092 }
00093 
00094 struct refSelStr {
00095     std::string segid;
00096     int resid;
00097     
00098     refSelStr() {} ;
00099     refSelStr(std::string newSegID, int newResid) {
00100         segid = newSegID;
00101         resid = newResid;
00102     }
00103 } ;
00104 
00105 typedef std::vector<refSelStr> refSelStrVec ;
00106 typedef std::map<Real, refSelStrVec> refSelStrMap ;
00107 typedef std::pair<Real,refSelStrVec> refSelStrPair ;
00108 
00109 void Molecule::prepare_qm(const char *pdbFileName, 
00110                                Parameters *params, ConfigList *cfgList) {
00111 #ifdef MEM_OPT_VERSION
00112     NAMD_die("QMMM interface is not supported in memory optimized builds.");
00113 #else
00114     
00115     PDB *pdbP;      //  Pointer to PDB object to use
00116     
00117     qmNumQMAtoms = 0;
00118     
00119     qmNumGrps = 0 ;
00120     
00121     iout << iINFO << "Using the following PDB file for QM parameters: " << 
00122     pdbFileName << "\n" << endi;
00123     if (pdbFileName)
00124         pdbP = new PDB(pdbFileName);
00125     else
00126         iout << "pdbFile name not defined!\n\n" << endi ;
00127     
00128     if (pdbP->num_atoms() != numAtoms)
00129     {
00130        NAMD_die("Number of atoms in QM parameter PDB file doesn't match coordinate PDB");
00131     }
00132     
00133     qmElementArray = new char*[numAtoms] ;
00134     
00135     qmAtomGroup = new Real[numAtoms] ;
00136     
00137     BigReal bondVal;
00138     Real atmGrp;
00139     
00140     std::set<Real> atmGrpSet;
00141     
00142     std::vector<Real> grpChrgVec;
00143     
00144     // Value in the bond column fro QM-MM bonds.
00145     std::vector<BigReal> qmBondValVec;
00146     // Identifiers of different QM regions
00147     std::vector<Real> qmGroupIDsVec;
00148     // This maps the qm group ID with the serail index of this group in 
00149     // various arrays.
00150     std::map<Real,int> qmGrpIDMap ;
00151     
00152     std::vector<int> qmAtmIndxVec, qmGrpSizeVec;
00153     
00154     // Used to parse user selection in different places.
00155     std::vector<std::string> strVec;
00156     
00157     qmNoPC = simParams->qmNoPC;
00158     
00159     // We set to zero by default, So there is no need for extra processing.
00160     qmPCFreq = 0;
00161     if (simParams->qmPCSelFreq > 1)
00162         qmPCFreq = simParams->qmPCSelFreq;
00163     
00164     
00167     
00168     
00169     qmLSSTotalNumAtms = 0;
00170     SortedArray< qmSolvData> lssGrpRes;
00171     std::map<Real,std::vector<int> > lssGrpRefIDs;
00172     refSelStrMap lssRefUsrSel;
00173     int totalNumRefAtms = 0;
00174     int lssClassicResIndx = 0 ;
00175     int lssCurrClassResID = -1 ;
00176     char lssCurrClassResSegID[5];
00177     if (simParams->qmLSSOn) {
00178         DebugM(4, "LSS is ON! Processing QM solvent.\n") ;
00179         
00180         if (resLookup == NULL)
00181             NAMD_die("We need residue data to conduct LSS.");
00182          
00183         if (simParams->qmLSSMode == QMLSSMODECOM) {
00184             
00185             StringList *current = cfgList->find("QMLSSRef");
00186             for ( ; current; current = current->next ) {
00187                 
00188                 strVec = split( std::string(current->data) , " ");
00189                 
00190                 if (strVec.size() != 3 ) {
00191                     iout << iERROR << "Format error in QM LSS size: " 
00192                     << current->data
00193                     << "\n" << endi;
00194                     NAMD_die("Error processing QM information.");
00195                 }
00196                 
00197                 std::stringstream storConv ;
00198                 
00199                 storConv << strVec[0] ;
00200                 Real grpID ;
00201                 storConv >> grpID;
00202                 if (storConv.fail()) {
00203                     iout << iERROR << "Error parsing QMLSSRef selection: " 
00204                     << current->data
00205                     << "\n" << endi;
00206                     NAMD_die("Error processing QM information.");
00207                 }
00208                 
00209                 std::string segName = strVec[1].substr(0,4);
00210                 
00211                 storConv.clear() ;
00212                 storConv << strVec[2];
00213                 int resID ;
00214                 storConv >> resID;
00215                 if (storConv.fail()) {
00216                     iout << iERROR << "Error parsing QMLSSRef selection: " 
00217                     << current->data
00218                     << "\n" << endi;
00219                     NAMD_die("Error processing QM information.");
00220                 }
00221                 
00222                 auto it = lssRefUsrSel.find(grpID) ;
00223                 if (it == lssRefUsrSel.end())
00224                     lssRefUsrSel.insert(refSelStrPair(grpID,refSelStrVec()));
00225                 
00226                 lssRefUsrSel[grpID].push_back(refSelStr(segName, resID));
00227             }
00228             
00229             for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
00230                 iout << iINFO << "LSS user selection for COM composition of group "
00231                 << it->first << ":\n" << endi ;
00232                 for (int i=0; i<it->second.size();i++) {
00233                     iout << iINFO << "Segment " << it->second[i].segid 
00234                     << " ; residue " << it->second[i].resid << "\n" << endi ;
00235                 }
00236             }
00237         }
00238     }
00239     
00240     
00241     
00244     
00245     
00246     for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
00247         
00248         // If we are looking for QM-MM bonds, then we need to store extra info.
00249         if (simParams->qmBondOn) {
00250             
00251             // We store both the qm group and the bond value
00252             if ( strcmp(simParams->qmColumn,"beta") == 0 ){
00253                 atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
00254                 
00255                 bondVal = pdbP->atom(atmInd)->occupancy() ;
00256             }
00257             else {
00258                 atmGrp = pdbP->atom(atmInd)->occupancy() ;
00259                 
00260                 bondVal = pdbP->atom(atmInd)->temperaturefactor() ;
00261             }
00262             
00263             qmBondValVec.push_back(bondVal);
00264         }
00265         else {
00266             
00267             if ( strcmp(simParams->qmColumn,"beta") == 0 ){
00268                 atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
00269             }
00270             else {
00271                 atmGrp = pdbP->atom(atmInd)->occupancy() ;
00272             }
00273         }
00274         
00275         // We store all atom QM Group IDs in an array for 
00276         // future transport to all PEs.
00277         qmAtomGroup[atmInd] = atmGrp;
00278         
00279         // We store all atom's elements for quick access in the QM code.
00280         // Elements are used to tell the QM code the atom type.
00281         qmElementArray[atmInd] = new char[3];
00282         strncpy(qmElementArray[atmInd],pdbP->atom(atmInd)->element(),3);
00283         
00284         // For QM atoms, we keep global atom indices and parse different QM Group
00285         // IDs, keeping track of each groups size
00286         if (atmGrp > 0){
00287             
00288             if (simParams->fixedAtomsOn) {
00289                 if ( fixedAtomFlags[atmInd] == 1 ) {
00290                     iout << iERROR << "QM atom cannot be fixed in space!\n" 
00291                     << endi;
00292                     NAMD_die("Error processing QM information.");
00293                 }
00294             }
00295             
00296             // Changes the VdW type of QM atoms.
00297             // This is may be used to counter balance the overpolarization 
00298             // that QM atoms sufer.
00299             if (simParams->qmVDW) {
00300                 // Creating a new type string
00301                 std::string qmType(qmElementArray[atmInd]) ;
00302                 // Erases empty spaces
00303                 qmType.erase(
00304                         std::remove_if(qmType.begin(), qmType.end(), isspace ),
00305                     qmType.end());
00306                 // pre-pends a "q" to the element name.
00307                 qmType = std::string("q") + qmType;
00308             
00309 //                 iout << "QM VdW type: " << atoms[atmInd].vdw_type 
00310 //                 << " atom type: " << atomNames[atmInd].atomtype 
00311 //                 << " new type "<< qmType << "\n" << endi;
00312                 
00313                 /*  Look up the vdw constants for this atom    */
00314                 // I am casting a non-const here because the function doesn't actually
00315                 // change the string value, but it doesn't ask for a const char* as 
00316                 // an argument.
00317                 params->assign_vdw_index(const_cast<char*>( qmType.c_str()), 
00318                                          &(atoms[atmInd]));
00319                 
00320 //                 iout << "----> new VdW type: " << atoms[atmInd].vdw_type << "\n" << endi;
00321             }
00322             
00323             // Indexes all global indices of QM atoms.
00324             qmAtmIndxVec.push_back(atmInd);
00325             
00326             auto retVal = atmGrpSet.insert(atmGrp);
00327             
00328             // If a new QM group is found, store the reverse reference in a map
00329             // and increase the total count.
00330             if (retVal.second) {
00331                 // This mak makes the reverse identification from the group ID
00332                 // to the sequential order in which each group was first found.
00333                 qmGrpIDMap.insert(std::pair<BigReal,int>(atmGrp, qmNumGrps));
00334                 
00335                 // This vector keeps track of the group ID for each group
00336                 qmGroupIDsVec.push_back(atmGrp);
00337                 
00338                 // This counter is used to keep track of the sequential order in
00339                 // which QM groups were first seen in the reference PDB file.
00340                 qmNumGrps++ ;
00341                 
00342                 // If this is a new QM group, initialize a new variable in the 
00343                 // vector to keep track of the number of atoms in this group.
00344                 qmGrpSizeVec.push_back(1);
00345                 
00346                 // For each new QM group, a new entry in the total charge
00347                 // vector is created
00348                 grpChrgVec.push_back( atoms[atmInd].charge );
00349             }
00350             else {
00351                 // If the QM group has already been seen, increase the count
00352                 // of atoms in that group.
00353                 qmGrpSizeVec[ qmGrpIDMap[atmGrp] ]++ ;
00354                 
00355                 // If the QM group exists, adds current atom charge to total charge.
00356                 grpChrgVec[ qmGrpIDMap[atmGrp] ] += atoms[atmInd].charge;
00357             }
00358             
00359             // In case we are using live solvent selection
00360             if(simParams->qmLSSOn) {
00361                 
00362                 int resID = pdbP->atom(atmInd)->residueseq();
00363                 char segName[5];
00364                 strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
00365                 
00366                 int resSize = get_residue_size(segName,resID);
00367                 
00368                 int i =-1, end =-1;
00369                 
00370                 resLookup->lookup(segName,resID,&i, &end);
00371                 
00372                 if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
00373                            simParams->qmLSSResname, 4) == 0) {
00374                     
00375                     // We try to insert every residue from every atom
00376                     qmSolvData *resP = lssGrpRes.find(qmSolvData(resID, i, resSize));
00377                     
00378                     if (resP != NULL) {
00379                         resP->atmIDs.push_back(atmInd) ;
00380 //                         DebugM(3, "Existing residue " << resP->resID 
00381 //                         << " from segID " << resP->segName
00382 //                         << " received atom "
00383 //                         << atmInd << "\n" );
00384                     }
00385                     else {
00386                         lssGrpRes.insert(qmSolvData(resID, i, resSize, segName, atmGrp));
00387 //                         DebugM(3, lssGrpRes.size() << ") Inserted new residue: " 
00388 //                         << resID << " from segID " << segName
00389 //                         << " with atom "
00390 //                         << i << "\n" ) ;
00391                     }
00392                 }
00393                 else {
00394                     // We store the QM atoms, per group, which are NOT part of
00395                     // solvent molecules.
00396                     
00397                     // Checks if we have a vector for this QM group.
00398                     auto itNewGrp = lssGrpRefIDs.find(atmGrp) ;
00399                     if (itNewGrp == lssGrpRefIDs.end()) {
00400                         lssGrpRefIDs.insert(std::pair<Real, std::vector<int> >(
00401                             atmGrp, std::vector<int>() ) );
00402                     }
00403                     
00404                     switch (simParams->qmLSSMode)
00405                     {
00406                     
00407                     case QMLSSMODECOM:
00408                     {
00409                         // If we are using COM selection, checks if the atom
00410                         // is part of the user selected residues
00411                         for(int i=0; i<lssRefUsrSel[atmGrp].size(); i++) {
00412                             if (lssRefUsrSel[atmGrp][i].resid == resID &&
00413                                 strncmp(lssRefUsrSel[atmGrp][i].segid.c_str(),
00414                                         segName,5) == 0 ) {
00415                                 
00416                                 lssGrpRefIDs[atmGrp].push_back(atmInd);
00417                                 totalNumRefAtms++;
00418                             }
00419                         }
00420                         
00421                     } break;
00422                     
00423                     case QMLSSMODEDIST:
00424                     {
00425                     
00426                         lssGrpRefIDs[atmGrp].push_back(atmInd);
00427                         totalNumRefAtms++;
00428                     
00429                     } break;
00430                     
00431                     }
00432                 }
00433                 
00434             }
00435             
00436         }
00437         else if (atmGrp == 0) {
00438             
00439             if(simParams->qmLSSOn) {
00440                 
00441                 if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
00442                            simParams->qmLSSResname, 4) == 0) {
00443                     
00444                     int resID = pdbP->atom(atmInd)->residueseq();
00445                     char segName[5];
00446                     strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
00447                     
00448                     if (lssCurrClassResID < 0) {
00449                         lssCurrClassResID = resID ;
00450                         strncpy(lssCurrClassResSegID, segName,5);
00451                         lssClassicResIndx = 0;
00452                     }
00453                     else if (lssCurrClassResID != resID ||
00454                             strcmp(lssCurrClassResSegID, segName) != 0 ) {
00455                         lssCurrClassResID = resID ;
00456                         strncpy(lssCurrClassResSegID, segName,5);
00457                         lssClassicResIndx++;
00458                     }
00459                     
00460                     qmClassicSolv.insert(std::pair<int,int>(atmInd,lssClassicResIndx));
00461                     
00462                 }
00463             }
00464         }
00465         else if(atmGrp < 0) {
00466             iout << iERROR << "QM group ID cannot be less than zero!\n" << endi;
00467             NAMD_die("Error processing QM information.");
00468         }
00469     }
00470     
00471     // Sanity check
00472     if (simParams->qmLSSOn) {
00473         if (lssGrpRes.size() == 0)
00474             NAMD_die("LSS was activated but there are no QM solvent molecules!\n");
00475         
00476         for (auto it=lssRefUsrSel.begin(); it!=lssRefUsrSel.end(); it++) {
00477             auto itTarget = qmGrpIDMap.find(it->first);
00478             if (itTarget == qmGrpIDMap.end()) {
00479                 iout << iERROR << "QM group ID for LSS could not be found in input!"
00480                 << " QM group ID: " << it->first << "\n" << endi;
00481                 NAMD_die("Error processing QM information.");
00482             }
00483         }
00484         
00485         DebugM(3,"We found " << lssClassicResIndx << " classical solvent molecules totalling "
00486         << qmClassicSolv.size() << " atoms.\n" );
00487         
00488     }
00489     
00490     qmNumQMAtoms = qmAtmIndxVec.size();
00491     
00492     if (qmNumQMAtoms == 0)
00493         NAMD_die("No QM atoms were found in this QM simulation!") ;
00494     
00495     iout << iINFO << "Number of QM atoms (excluding Dummy atoms): " << 
00496     qmNumQMAtoms << "\n" << endi;
00497     
00498     qmAtmChrg = new Real[qmNumQMAtoms] ;
00499     qmAtmIndx = new int[qmNumQMAtoms] ;
00500     for (int i=0; i<qmNumQMAtoms; i++) {
00501         // qmAtmChrg gets initialized with the PSF charges at the end of this
00502         // function, but values may change as QM steps are taken.
00503         qmAtmChrg[i] = 0;  
00504         qmAtmIndx[i] = qmAtmIndxVec[i] ;
00505     }
00506     
00507     // This map relates the QM group index with a vector of pairs
00508     // of bonded MM-QM atoms (with the bonded atoms ins this order: 
00509     // MM first, QM second).
00510     std::map<int, std::vector<std::pair<int,int> > > qmGrpIDBonds ;
00511     int bondCounter = 0;
00512     if (simParams->qmBondOn) {
00513         
00514         // Checks all bonds for QM-MM pairs.
00515         // Makes sanity checks against QM-QM pairs and MM-MM pairs that
00516         // are flagged by the user to be bonds between QM and MM regions.
00517         // QM-QM bonds will be removed in another function.
00518         for (int bndIt = 0; bndIt < numBonds; bndIt++) {
00519             
00520             bond curr = bonds[bndIt] ;
00521             
00522             // In case either atom has a non-zero
00523             if ( qmBondValVec[curr.atom1] != 0 &&
00524                  qmBondValVec[curr.atom2] != 0 ) {
00525                 
00526                 // Sanity checks (1 of 2)
00527                 if (qmAtomGroup[curr.atom1] != 0 &&
00528                     qmAtomGroup[curr.atom2] != 0) {
00529                     iout << iERROR << "Atoms " << curr.atom1 << " and " << 
00530                     curr.atom2 << " are assigned as QM atoms.\n" << endi;
00531                     NAMD_die("Error in QM-MM bond assignment.") ;
00532                 }
00533                 
00534                 // Sanity checks (2 of 2)
00535                 if (qmAtomGroup[curr.atom1] == 0 &&
00536                     qmAtomGroup[curr.atom2] == 0) {
00537                     iout << iERROR << "Atoms " << curr.atom1 << " and " << 
00538                     curr.atom2 << " are assigned as MM atoms.\n" << endi;
00539                     NAMD_die("Error in QM-MM bond assignment.") ;
00540                 }
00541                 
00542                 int currGrpID = 0;
00543                 std::pair<int,int> newPair(0,0);
00544                 
00545                 // We create a QM-MM pair with the MM atom first
00546                 if (qmAtomGroup[curr.atom1] != 0) {
00547                     newPair = std::pair<int,int>(curr.atom2, curr.atom1) ;
00548                     currGrpID = qmAtomGroup[curr.atom1];
00549                 } else {
00550                     newPair = std::pair<int,int>(curr.atom1, curr.atom2) ;
00551                     currGrpID = qmAtomGroup[curr.atom2];
00552                 } 
00553                 
00554                 int grpIndx = qmGrpIDMap[currGrpID] ;
00555                 
00556                 // Stores bonds in vectors key-ed by the QM group ID of the QM atom.
00557                 auto retIter = qmGrpIDBonds.find(grpIndx) ;
00558                 // In case thi QM-MM bonds belong to a QM group we have not seen 
00559                 // yet, stores a new vector in the map.
00560                 if (retIter == qmGrpIDBonds.end()) {
00561                     qmGrpIDBonds.insert(std::pair<BigReal,std::vector<std::pair<int,int> > >(
00562                     grpIndx, std::vector<std::pair<int,int> >() ) );
00563                 }
00564                 
00565                 qmGrpIDBonds[grpIndx].push_back( newPair );
00566                 
00567                 bondCounter++ ;
00568             }
00569             
00570             
00571         }
00572         
00573 //         iout << "Finished processing QM-MM bonds.\n" << endi ;
00574         
00575         if (bondCounter == 0)
00576             iout << iWARN << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
00577         else
00578             iout << iINFO << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
00579         
00580     }
00581     
00582     // Initializes several arrays used to setup the QM simulation.
00583     qmNumBonds = bondCounter ;
00584     
00585     qmMMBond = new int*[qmNumBonds];
00586     
00587     qmDummyBondVal = new BigReal[qmNumBonds];
00588     
00589     qmMMChargeTarget = new int*[qmNumBonds] ;
00590     qmMMNumTargs = new int[qmNumBonds] ;
00591     
00592     qmDummyElement = new char*[qmNumBonds] ;
00593     
00594     
00595     qmNumGrps = atmGrpSet.size();
00596     
00597     qmGrpSizes = new int[qmNumGrps];
00598     
00599     qmGrpID = new Real[qmNumGrps];
00600     
00601     qmGrpChrg = new Real[qmNumGrps];
00602     
00603     qmGrpMult = new Real[qmNumGrps] ;
00604     
00605     qmGrpNumBonds = new int[qmNumGrps];
00606     
00607     qmGrpBonds = new int*[qmNumGrps];
00608     qmMMBondedIndx = new int*[qmNumGrps];
00609     
00610     
00613     
00614     
00615     // We first initialize the multiplicity vector with 
00616     // default values, then populate it with user defined values.
00617     for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
00618         qmGrpMult[grpIter] = 0;
00619     }
00620     
00621     int multCount = 0;
00622     StringList *current = cfgList->find("QMMult");
00623     for ( ; current; current = current->next ) {
00624         
00625         auto strVec = split( std::string(current->data) , " ");
00626         
00627         if (strVec.size() != 2 ) {
00628             iout << iERROR << "Format error in QM Multiplicity string: " 
00629             << current->data
00630             << "\n" << endi;
00631             NAMD_die("Error processing QM information.");
00632         }
00633         
00634         std::stringstream storConv ;
00635         
00636         storConv << strVec[0] ;
00637         Real grpID ;
00638         storConv >> grpID;
00639         if (storConv.fail()) {
00640             iout << iERROR << "Error parsing QMMult selection: " 
00641             << current->data
00642             << "\n" << endi;
00643             NAMD_die("Error processing QM information.");
00644         }
00645         
00646         storConv.clear() ;
00647         storConv << strVec[1];
00648         Real multiplicity ;
00649         storConv >> multiplicity;
00650         if (storConv.fail()) {
00651             iout << iERROR << "Error parsing QMMult selection: " 
00652             << current->data
00653             << "\n" << endi;
00654             NAMD_die("Error processing QM information.");
00655         }
00656         
00657         auto it = qmGrpIDMap.find(grpID);
00658         
00659         if (it == qmGrpIDMap.end()) {
00660             iout << iERROR << "Error parsing QMMult selection. Could not find QM group ID: " 
00661             << grpID
00662             << "\n" << endi;
00663             NAMD_die("Error processing QM information.");
00664         }
00665         else {
00666             iout << iINFO << "Applying user defined multiplicity "
00667             << multiplicity << " to QM group ID " << grpID << "\n" << endi;
00668             qmGrpMult[it->second] = multiplicity;
00669         }
00670         
00671         multCount++;
00672     }
00673     
00674     if (multCount != qmNumGrps && simParams->qmFormat == QMFormatORCA) {
00675         NAMD_die("ORCA neds multiplicity values for all QM regions.");
00676     }
00677     
00678     
00681     
00682     
00683     for (size_t grpIndx = 0; grpIndx < qmGrpSizeVec.size(); grpIndx++) {
00684         
00685         bool nonInteger = true;
00686         if ((fabsf(roundf(grpChrgVec[grpIndx]) - grpChrgVec[grpIndx]) <= 0.001f) ) {
00687             grpChrgVec[grpIndx] = roundf(grpChrgVec[grpIndx]) ;
00688             nonInteger = false;
00689         }
00690         
00691         iout << iINFO << grpIndx + 1 << ") Group ID: " << qmGroupIDsVec[grpIndx]
00692         << " ; Group size: " << qmGrpSizeVec[grpIndx] << " atoms"
00693         << " ; Total PSF charge: " << grpChrgVec[grpIndx] << "\n" << endi ;
00694         
00695         if (nonInteger && simParams->PMEOn)
00696             NAMD_die("QM atoms do not add up to a whole charge, which is needed for PME.") ;
00697     }
00698     
00699     int chrgCount = 0;
00700     current = cfgList->find("QMCharge");
00701     for ( ; current; current = current->next ) {
00702         
00703         auto strVec = split( std::string(current->data) , " ");
00704         
00705         if (strVec.size() != 2 ) {
00706             iout << iERROR << "Format error in QM Charge string: " 
00707             << current->data
00708             << "\n" << endi;
00709             NAMD_die("Error processing QM information.");
00710         }
00711         
00712         std::stringstream storConv ;
00713         
00714         storConv << strVec[0] ;
00715         Real grpID ;
00716         storConv >> grpID;
00717         if (storConv.fail()) {
00718             iout << iERROR << "Error parsing QMCharge selection: " 
00719             << current->data
00720             << "\n" << endi;
00721             NAMD_die("Error processing QM information.");
00722         }
00723         
00724         storConv.clear() ;
00725         storConv << strVec[1];
00726         Real charge ;
00727         storConv >> charge;
00728         if (storConv.fail()) {
00729             iout << iERROR << "Error parsing QMCharge selection: " 
00730             << current->data
00731             << "\n" << endi;
00732             NAMD_die("Error processing QM information.");
00733         }
00734         
00735         auto it = qmGrpIDMap.find(grpID);
00736         
00737         if (it == qmGrpIDMap.end()) {
00738             iout << iERROR << "Error parsing QMCharge selection. Could not find QM group ID: " 
00739             << grpID
00740             << "\n" << endi;
00741             NAMD_die("Error processing QM information.");
00742         }
00743         else {
00744             iout << iINFO << "Found user defined charge "
00745             << charge << " for QM group ID " << grpID << ". Will ignore PSF charge.\n" << endi;
00746             grpChrgVec[it->second] = charge;
00747         }
00748         
00749         chrgCount++;
00750     }
00751     
00752     simParams->qmMOPACAddConfigChrg = false;
00753     // Checks if QM group charges can be defined for MOPAC.
00754     // Since charge can be defined in QMConfigLine, we need extra logic.
00755     if (simParams->qmFormat == QMFormatMOPAC) {
00756         
00757         // Checks if group charge was supplied in config line for MOPAC.
00758         std::string::size_type chrgLoc = std::string::npos ;
00759         
00760         // This will hold a sting with the first user supplied configuration line for MOPAC.
00761         std::string configLine(cfgList->find("QMConfigLine")->data) ;
00762         chrgLoc = configLine.find("CHARGE") ;
00763         
00764         if ( chrgLoc != std::string::npos ) {
00765             iout << iINFO << "Found user defined charge in command line. This \
00766 will be used for all QM regions and will take precedence over all other charge \
00767 definitions.\n" << endi;
00768         }
00769         else if ( chrgLoc == std::string::npos && (simParams->qmChrgFromPSF || chrgCount == qmNumGrps) ) {
00770         // If no charge was defined in the configuration line, gets from PSF and/or 
00771         // from user defined charges (through the QMCharge keyword).
00772             simParams->qmMOPACAddConfigChrg = true;
00773         }
00774         else
00775         {
00776         // If we could nont find a charge definition in the config line AND
00777         // no specific charge was selected for each QM region through QMCharge AND
00778         // the QMChargeFromPSF was not turned ON, the we stop NAMD and scream at the user.
00779 //         if ( chrgLoc == std::string::npos && (chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) 
00780             iout << iERROR << "Could not find charge for all QM groups. For MOPAC, \
00781 charges can be defined through QMConfigLine, QMCharge or QMChargeFromPSF keywords.\n" << endi;
00782             NAMD_die("Error processing QM information.");
00783         }
00784         
00785     }
00786     
00787     if (simParams->qmFormat == QMFormatORCA ) {
00788         if ((chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) {
00789             // If we are not supposed to get charges from the PSF, and not 
00790             // enough charges were set to cover all QM regions, 
00791             // we stop NAMD and scream at the user.
00792             iout << iERROR << "Could not find charge for all QM groups. For ORCA, \
00793 charges can be defined through QMCharge or QMChargeFromPSF keywords.\n" << endi;
00794             NAMD_die("Error processing QM information.");
00795         }
00796     }
00797     
00798     // If mechanichal embedding was requested but we have QM-MM bonds, we need 
00799     // to send extra info to ComputeQM to preserve calculation speed.
00800     if (qmNumBonds > 0 && qmNoPC) {
00801         qmMeNumBonds = qmNumBonds;
00802         qmMeMMindx = new int[qmMeNumBonds] ;
00803         qmMeQMGrp = new Real[qmMeNumBonds] ;
00804     }
00805     else {
00806         qmMeNumBonds = 0 ;
00807         qmMeMMindx = NULL;
00808         qmMeQMGrp = NULL;
00809     }
00810     
00811     
00814     
00815     
00816     bondCounter = 0;
00817     for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
00818         
00819         qmGrpID[grpIter] = qmGroupIDsVec[grpIter] ;
00820         
00821         qmGrpSizes[grpIter] = qmGrpSizeVec[grpIter] ;
00822         
00823         qmGrpChrg[grpIter] = grpChrgVec[grpIter];
00824         
00825 //         iout << "Loaded " << qmGrpSizes[grpIter] << " atoms to this group.\n" << endi;
00826         
00827         int currNumbBonds = qmGrpIDBonds[grpIter].size() ;
00828         
00829         // Assigns the number of bonds that the current QM group has.
00830         qmGrpNumBonds[grpIter] = currNumbBonds;
00831         
00832         if (currNumbBonds > 0) {
00833             
00834             qmGrpBonds[grpIter] = new int[currNumbBonds];
00835             qmMMBondedIndx[grpIter] = new int[currNumbBonds];
00836             
00837             for (int bndIter = 0; bndIter<currNumbBonds; bndIter++) {
00838                 
00839                 // Adds the bonds to the overall sequential list.
00840                 qmMMBond[bondCounter] = new int[2] ;
00841                 qmMMBond[bondCounter][0] = qmGrpIDBonds[grpIter][bndIter].first ;
00842                 qmMMBond[bondCounter][1] = qmGrpIDBonds[grpIter][bndIter].second ;
00843                 
00844                 // For the current QM group, and the current bond, gets the bond index.
00845                 qmGrpBonds[grpIter][bndIter] =  bondCounter;
00846                 
00847                 // For the current QM group, and the current bond, gets the MM atom.
00848                 qmMMBondedIndx[grpIter][bndIter] = qmGrpIDBonds[grpIter][bndIter].first ;
00849                 
00850                 // Assign the default value of dummy element
00851                 qmDummyElement[bondCounter] = new char[3];
00852                 strcpy(qmDummyElement[bondCounter],"H\0");
00853                 
00854                 // Determines the distance that will separate the new Dummy atom
00855                 // and the Qm atom to which it will be bound.
00856                 bondVal = 0;
00857                 if (simParams->qmBondDist) {
00858                     if (qmBondValVec[qmMMBond[bondCounter][0]] != 
00859                         qmBondValVec[qmMMBond[bondCounter][1]]
00860                     ) {
00861                         iout << iERROR << "qmBondDist is ON but the values in the bond column are different for atom "
00862                         << qmMMBond[bondCounter][0] << " and " << qmMMBond[bondCounter][1]
00863                         << ".\n" << endi ;
00864                         NAMD_die("Error in QM-MM bond processing.");
00865                     }
00866                     
00867                     bondVal = qmBondValVec[qmMMBond[bondCounter][0]] ;
00868                 } 
00869                 else {
00870                     
00871                     if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"C") == 0 ) {
00872                         bondVal = 1.09 ;
00873                     }
00874                     else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"O") == 0 ) {
00875                         bondVal = 0.98 ;
00876                     }
00877                     else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"N") == 0 ) {
00878                         bondVal = 0.99 ;
00879                     }
00880                     else {
00881                         iout << iERROR << "qmBondDist is OFF but the QM atom has no default distance value. Atom "
00882                         << qmMMBond[bondCounter][1] << ", with element: " 
00883                         << qmElementArray[qmMMBond[bondCounter][1]] 
00884                         <<  ".\n" << endi ;
00885                         NAMD_die("Error in QM-MM bond processing.");
00886                     }
00887                     
00888                 }
00889                 
00890                 iout << iINFO << "MM-QM pair: " << qmMMBond[bondCounter][0] << ":"
00891                 << qmMMBond[bondCounter][1] 
00892                 << " -> Value (distance or ratio): " << bondVal
00893                 << " (QM Group " << grpIter << " ID " << qmGrpID[grpIter] << ")"
00894                 << "\n" << endi ;
00895                 
00896                 qmDummyBondVal[bondCounter] = bondVal;
00897                 
00898                 // In case we are preparing for a mechanical embedding simulation
00899                 // with no point charges, populate the following vectors
00900                 if (qmMeNumBonds > 0) {
00901                     qmMeMMindx[bondCounter] = qmMMBond[bondCounter][0];
00902                     qmMeQMGrp[bondCounter] = qmGrpID[grpIter];
00903                 }
00904                 
00905                 bondCounter++ ;
00906             }
00907         }
00908         
00909     }
00910     
00911     
00914     
00915     current = NULL;
00916     if (qmNumBonds > 0)
00917         current = cfgList->find("QMLinkElement");
00918     
00919     int numParsedLinkElem = 0;
00920     for ( ; current != NULL; current = current->next ) {
00921         
00922         DebugM(3,"Parsing link atom element: " << current->data << "\n" );
00923         
00924         strVec = split( std::string(current->data) , " ");
00925         
00926         // We need the two atoms that compose the QM-MM bonds and 
00927         // then the element.
00928         if (strVec.size() != 3 && qmNumBonds > 1) {
00929             iout << iERROR << "Format error in QM link atom element selection: " 
00930             << current->data
00931             << "\n" << endi;
00932             NAMD_die("Error processing QM information.");
00933         }
00934         
00935         // If there is only one QM-MM bond, we can accept the element only.
00936         if (strVec.size() != 1 && strVec.size() != 3 && qmNumBonds == 1) {
00937             iout << iERROR << "Format error in QM link atom element selection: " 
00938             << current->data
00939             << "\n" << endi;
00940             NAMD_die("Error processing QM information.");
00941         }
00942         
00943         std::stringstream storConv ;
00944         int bondAtom1, bondAtom2;
00945         std::string linkElement ;
00946         
00947         if (strVec.size() == 1) {
00948             linkElement = strVec[0].substr(0,2);
00949         }
00950         else if (strVec.size() == 3) {
00951             
00952             storConv << strVec[0] ;
00953             storConv >> bondAtom1;
00954             if (storConv.fail()) {
00955                 iout << iERROR << "Error parsing link atom element selection: " 
00956                 << current->data
00957                 << "\n" << endi;
00958                 NAMD_die("Error processing QM information.");
00959             }
00960             
00961             storConv.clear() ;
00962             storConv << strVec[1];
00963             storConv >> bondAtom2;
00964             if (storConv.fail()) {
00965                 iout << iERROR << "Error parsing link atom element selection: " 
00966                 << current->data
00967                 << "\n" << endi;
00968                 NAMD_die("Error processing QM information.");
00969             }
00970             
00971             linkElement = strVec[2].substr(0,2);
00972         }
00973         
00974         numParsedLinkElem++;
00975         
00976         if (numParsedLinkElem>qmNumBonds) {
00977             iout << iERROR << "More link atom elements were set than there"
00978             " are QM-MM bonds.\n" << endi;
00979             NAMD_die("Error processing QM information.");
00980         }
00981         
00982         int bondIter;
00983         
00984         if (strVec.size() == 1) {
00985             bondIter = 0;
00986         }
00987         else if (strVec.size() == 3) {
00988             
00989             Bool foundBond = false;
00990             
00991             for (bondIter=0; bondIter<qmNumBonds; bondIter++) {
00992                 
00993                 if ( (qmMMBond[bondIter][0] == bondAtom1 &&
00994                     qmMMBond[bondIter][1] == bondAtom2 ) ||
00995                     (qmMMBond[bondIter][0] == bondAtom2 &&
00996                     qmMMBond[bondIter][1] == bondAtom1 ) ) {
00997                     
00998                     foundBond = true;
00999                     break;
01000                 }
01001             }
01002             
01003             if (! foundBond) {
01004                 iout << iERROR << "Error parsing link atom element selection: " 
01005                 << current->data
01006                 << "\n" << endi;
01007                 iout << iERROR << "Atom pair was not found to be a QM-MM bond.\n"
01008                 << endi;
01009                 NAMD_die("Error processing QM information.");
01010             }
01011         }
01012         
01013         strcpy(qmDummyElement[bondIter],linkElement.c_str());
01014         qmDummyElement[bondIter][2] = '\0';
01015         
01016         iout << iINFO << "Applying user defined link atom element "
01017         << qmDummyElement[bondIter] << " to QM-MM bond "
01018         << bondIter << ": " << qmMMBond[bondIter][0]
01019         << " - " << qmMMBond[bondIter][1]
01020         << "\n" << endi;
01021     }
01022     
01023     
01024     
01027     
01028     
01029     int32 **bondsWithAtomLocal = NULL ;
01030     int32 *byAtomSizeLocal = NULL;
01031     ObjectArena <int32 >* tmpArenaLocal = NULL;
01032     if (qmNumBonds > 0) {
01033         
01034         bondsWithAtomLocal = new int32 *[numAtoms];
01035         byAtomSizeLocal = new int32[numAtoms];
01036         tmpArenaLocal = new ObjectArena<int32>;
01037         
01038         //  Build the bond lists
01039        for (int i=0; i<numAtoms; i++)
01040        {
01041          byAtomSizeLocal[i] = 0;
01042        }
01043        for (int i=0; i<numRealBonds; i++)
01044        {
01045          byAtomSizeLocal[bonds[i].atom1]++;
01046          byAtomSizeLocal[bonds[i].atom2]++;
01047        }
01048        for (int i=0; i<numAtoms; i++)
01049        {
01050          bondsWithAtomLocal[i] = tmpArenaLocal->getNewArray(byAtomSizeLocal[i]+1);
01051          bondsWithAtomLocal[i][byAtomSizeLocal[i]] = -1;
01052          byAtomSizeLocal[i] = 0;
01053        }
01054        for (int i=0; i<numRealBonds; i++)
01055        {
01056          int a1 = bonds[i].atom1;
01057          int a2 = bonds[i].atom2;
01058          bondsWithAtomLocal[a1][byAtomSizeLocal[a1]++] = i;
01059          bondsWithAtomLocal[a2][byAtomSizeLocal[a2]++] = i;
01060        }
01061     }
01062     
01063     // In this loops we try to find other bonds in which the MM atoms from
01064     // QM-MM bonds may be involved. The other MM atoms (which we will call M2 and M3)
01065     // will be involved in charge manipulation. See ComputeQM.C for details.
01066     for (int qmBndIt = 0; qmBndIt < qmNumBonds; qmBndIt++) {
01067         
01068         // The charge targets are accumulated in a temporary vector and then 
01069         // transfered to an array that will be transmited to the ComputeQMMgr object.
01070         std::vector<int> chrgTrgt ;
01071         
01072         int MM1 = qmMMBond[qmBndIt][0], MM2, MM2BondIndx, MM3, MM3BondIndx;
01073         
01074         switch (simParams->qmBondScheme) {
01075             
01076             case QMSCHEMERCD:
01077             
01078             case QMSCHEMECS:
01079             {
01080                 // Selects ALL MM2 atoms.
01081                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
01082                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
01083                     
01084                     // Checks which of the atoms in the bond structure is the
01085                     // MM2 atom.
01086                     if (bonds[MM2BondIndx].atom1 == MM1)
01087                         MM2 = bonds[MM2BondIndx].atom2;
01088                     else
01089                         MM2 = bonds[MM2BondIndx].atom1;
01090                     
01091                     // In case the bonded atom is a QM atom,
01092                     // skips the index.
01093                     if (qmAtomGroup[MM2] > 0)
01094                         continue;
01095                     
01096                     chrgTrgt.push_back(MM2);
01097                 }
01098                 
01099             } break;
01100             
01101             case QMSCHEMEZ3:
01102             {
01103                 // Selects all MM3 atoms
01104                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
01105                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
01106                     
01107                     // Checks which of the atoms in the bond structure is the
01108                     // MM2 atom.
01109                     if (bonds[MM2BondIndx].atom1 == MM1)
01110                         MM2 = bonds[MM2BondIndx].atom2;
01111                     else
01112                         MM2 = bonds[MM2BondIndx].atom1;
01113                     
01114                     // In case the bonded atom is a QM atom,
01115                     // skips the index.
01116                     if (qmAtomGroup[MM2] > 0)
01117                         continue;
01118                     
01119                     for (int i=0; i<byAtomSizeLocal[MM2]; i++) {
01120                         MM3BondIndx = bondsWithAtomLocal[MM2][i] ;
01121                         
01122                         // Checks which of the atoms in the bond structure is the
01123                         // MM3 atom.
01124                         if (bonds[MM3BondIndx].atom1 == MM2)
01125                             MM3 = bonds[MM3BondIndx].atom2;
01126                         else
01127                             MM3 = bonds[MM3BondIndx].atom1;
01128                         
01129                         // In case the bonded atom is a QM atom,
01130                         // skips the index.
01131                         // We also keep the search from going back to MM1.
01132                         if (qmAtomGroup[MM3] > 0 || MM3 == MM1)
01133                             continue;
01134                         
01135                         chrgTrgt.push_back(MM3);
01136                     }
01137                     
01138                 }
01139                 
01140             };
01141             
01142             case QMSCHEMEZ2:
01143             {
01144                 // Selects all MM2 atoms
01145                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
01146                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
01147                     
01148                     // Checks which of the atoms in the bond structure is the
01149                     // MM2 atom.
01150                     if (bonds[MM2BondIndx].atom1 == MM1)
01151                         MM2 = bonds[MM2BondIndx].atom2;
01152                     else
01153                         MM2 = bonds[MM2BondIndx].atom1;
01154                     
01155                     // In case the bonded atom is a QM atom,
01156                     // skips the index.
01157                     if (qmAtomGroup[MM2] > 0)
01158                         continue;
01159                     
01160                     chrgTrgt.push_back(MM2);
01161                 }
01162                 
01163             };
01164             
01165             case QMSCHEMEZ1:
01166             {
01167                 // Selects all MM1 atoms
01168                 chrgTrgt.push_back(MM1);
01169             } break;
01170         }
01171         
01172         
01173         qmMMChargeTarget[qmBndIt] = new int[chrgTrgt.size()] ;
01174         qmMMNumTargs[qmBndIt] =  chrgTrgt.size();
01175         
01176         DebugM(3, "MM-QM bond " << qmBndIt << "; MM atom " 
01177         << qmMMBond[qmBndIt][0] << " conections: \n" );
01178         
01179         for (size_t i=0; i < chrgTrgt.size(); i++) {
01180             qmMMChargeTarget[qmBndIt][i] = chrgTrgt[i];
01181             DebugM(3,"MM Bonded to: " << chrgTrgt[i] << "\n" );
01182         }
01183         
01184         chrgTrgt.clear();
01185     }
01186     
01187     if (bondsWithAtomLocal != NULL)
01188         delete [] bondsWithAtomLocal;  bondsWithAtomLocal = 0;
01189     if (byAtomSizeLocal != NULL)
01190         delete [] byAtomSizeLocal;  byAtomSizeLocal = 0;
01191     if (tmpArenaLocal != NULL)
01192         delete tmpArenaLocal;  tmpArenaLocal = 0;
01193     
01194     
01197     
01198     
01199     if(simParams->qmLSSOn) {
01200         
01201         std::map<Real,int> grpLSSSize ;
01202         std::map<Real,int>::iterator itGrpSize;
01203         
01204         qmLSSTotalNumAtms = 0;
01205         qmLSSResidueSize = 0;
01206         
01207         if (simParams->qmLSSFreq == 0)
01208             qmLSSFreq = simParams->stepsPerCycle ;
01209         else
01210             qmLSSFreq = simParams->qmLSSFreq;
01211             
01212         #ifdef DEBUG_QM
01213         int resSize = -1;
01214         #endif
01215         
01216         std::map<Real, int> grpNumLSSRes;
01217         std::map<Real, int>::iterator itGrpNumRes;
01218         
01219         for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
01220             
01221             if (it->atmIDs.size() != it->size) {
01222                 iout << iERROR << "The number of atoms loaded for residue "
01223                 << it->resID << " does not match the expected for this residue type.\n" 
01224                 << endi;
01225                 NAMD_die("Error parsing data for LSS.");
01226             }
01227             
01228             qmLSSTotalNumAtms += it->size;
01229             
01230             #ifdef DEBUG_QM
01231             if (resSize < 0) resSize = it->size ;
01232             if (resSize > 0 and resSize != it->size) {
01233                 iout << iERROR << "The number of atoms loaded for residue "
01234                 << it->resID << " does not match previously loaded residues.\n" 
01235                 << endi;
01236                 NAMD_die("Error parsing data for LSS.");
01237             }
01238                 
01239 //             DebugM(3,"Residue " << it->resID << ": " << it->segName
01240 //             << " - from " << it->begAtmID << " with size "
01241 //             << it->size << " (QM ID: " << it->qmGrpID
01242 //             << ") has " << it->atmIDs.size() << " atoms: \n" ) ;
01243 //             for (int i=0; i<it->atmIDs.size(); i++) 
01244 //                 DebugM(3, it->atmIDs[i] << "\n" );
01245             #endif
01246             
01247             // Calculating total number of atoms per group
01248             itGrpSize = grpLSSSize.find(it->qmGrpID) ;
01249             if (itGrpSize != grpLSSSize.end())
01250                 itGrpSize->second += it->size;
01251             else
01252                 grpLSSSize.insert(std::pair<Real,int>(it->qmGrpID, it->size));
01253             
01254             // Calculating total number of solvent residues per group
01255             itGrpNumRes = grpNumLSSRes.find(it->qmGrpID) ;
01256             if (itGrpNumRes != grpNumLSSRes.end())
01257                 itGrpNumRes->second += 1;
01258             else
01259                 grpNumLSSRes.insert(std::pair<Real,int>(it->qmGrpID, 1));
01260         }
01261         
01262         qmLSSResidueSize = lssGrpRes.begin()->size ;
01263         
01264         qmLSSSize = new int[qmNumGrps];
01265         
01266         qmLSSIdxs = new int[qmLSSTotalNumAtms];
01267         int lssAtmIndx = 0;
01268         
01269         switch (simParams->qmLSSMode) {
01270         
01271         case QMLSSMODECOM:
01272         {
01273             
01274             qmLSSRefSize = new int[qmNumGrps];
01275             
01276             qmLSSMass = new Mass[qmLSSTotalNumAtms];
01277             
01278             qmLSSRefIDs = new int[totalNumRefAtms];
01279             qmLSSRefMass = new Mass[totalNumRefAtms];
01280             int lssRefIndx = 0;
01281             
01282             for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
01283                 
01284                 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
01285                 
01286                 if (itGrpSize != grpNumLSSRes.end())
01287                     qmLSSSize[grpIndx] =  itGrpSize->second;
01288                 else
01289                     qmLSSSize[grpIndx] =  0;
01290                 
01291                 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
01292                     
01293                     if (it->qmGrpID == qmGrpID[grpIndx]) {
01294                         for (int i=0; i<it->atmIDs.size(); i++) {
01295                             qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
01296                             qmLSSMass[lssAtmIndx] = atoms[it->atmIDs[i]].mass;
01297                             lssAtmIndx++;
01298                         }
01299                     }
01300                 }
01301                 
01302                 DebugM(4, "QM group " << qmGrpID[grpIndx]
01303                 << " has " << qmLSSSize[grpIndx] << " solvent molecules, "
01304                 << "totalling " << grpLSSSize[qmGrpID[grpIndx]]
01305                 << " atoms (check " << lssAtmIndx << ").\n" );
01306                 
01307                 qmLSSRefSize[grpIndx] = lssGrpRefIDs[qmGrpID[grpIndx]].size();
01308                 for(int i=0; i < qmLSSRefSize[grpIndx]; i++) {
01309                     qmLSSRefIDs[lssRefIndx] = lssGrpRefIDs[qmGrpID[grpIndx]][i];
01310                     qmLSSRefMass[lssRefIndx] =  atoms[qmLSSRefIDs[lssRefIndx]].mass;
01311                     lssRefIndx++;
01312                 }
01313                 
01314                 DebugM(4, "QM group " << qmGrpID[grpIndx]
01315                 << " has " << qmLSSRefSize[grpIndx] << " non-solvent atoms for LSS.\n" );
01316             }
01317             
01318         } break ;
01319         
01320         case QMLSSMODEDIST:
01321         {
01322             for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
01323                 
01324                 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
01325                 
01326                 if (itGrpSize != grpNumLSSRes.end())
01327                     qmLSSSize[grpIndx] =  itGrpSize->second;
01328                 else
01329                     qmLSSSize[grpIndx] =  0;
01330                 
01331                 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
01332                     
01333                     if (it->qmGrpID == qmGrpID[grpIndx]) {
01334                         for (int i=0; i<it->atmIDs.size(); i++) {
01335                             qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
01336                             lssAtmIndx++;
01337                         }
01338                     }
01339                 }
01340                 
01341                 DebugM(4, "QM group " << qmGrpID[grpIndx]
01342                 << " has " << qmLSSSize[grpIndx] << " solvent molecules, "
01343                 << "totalling " << grpLSSSize[qmGrpID[grpIndx]]
01344                 << " atoms (check " << lssAtmIndx << ").\n" );
01345             }
01346             
01347         } break ;
01348             
01349         }
01350     }
01351     
01352     
01355     
01356     
01357     PDB *customPCPDB;
01358     
01359     // In case we have a custom and fixed set of point charges for each QM group,
01360     // we process the files containing information.
01361     current = NULL;
01362     if (simParams->qmCustomPCSel) {
01363         current = cfgList->find("QMCustomPCFile");
01364     }
01365     
01366     std::map<Real,std::vector<int> > qmPCVecMap ;
01367     
01368     int numParsedPBDs = 0;
01369     for ( ; current != NULL; current = current->next ) {
01370         
01371         iout << iINFO << "Parsing QM Custom PC file " << current->data << "\n" << endi;
01372         customPCPDB = new PDB(current->data);
01373         
01374         if (customPCPDB->num_atoms() != numAtoms)
01375            NAMD_die("Number of atoms in QM Custom PC PDB file doesn't match coordinate PDB");
01376         
01377         std::vector< int > currPCSel ;
01378         Real currQMID = 0 ;
01379         int currGrpSize = 0 ;
01380         
01381         for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
01382             
01383             BigReal beta = customPCPDB->atom(atmInd)->temperaturefactor() ;
01384             BigReal occ = customPCPDB->atom(atmInd)->occupancy() ;
01385             
01386             if ( beta != 0 && occ != 0)
01387                 NAMD_die("An atom cannot be marked as QM and poitn charge simultaneously!");
01388             
01389             // If this is not a QM atom and 
01390             if (occ != 0) {
01391                 currPCSel.push_back(atmInd) ;
01392             }
01393             
01394             if (beta != 0) {
01395                 if (pdbP->atom(atmInd)->temperaturefactor() != beta)
01396                     NAMD_die("QM Group selection is different in reference PDB and Custom-PC PDB!");
01397                 
01398                 if (currQMID == 0) {
01399                     // If this is the first QM atom we find, keep the QM Group ID.
01400                     currQMID = beta;
01401                 }
01402                 else {
01403                     // For all other atoms, check if it is the same group. It must be!!
01404                     if (currQMID != beta)
01405                         NAMD_die("Found two different QM group IDs in this file!");
01406                 }
01407                 
01408                 currGrpSize++;
01409             }
01410             
01411         }
01412         
01413         if (currGrpSize != qmGrpSizeVec[ qmGrpIDMap[currQMID] ])
01414             NAMD_die("Reference PDB and Custom-PC PDB have different QM group sizes!") ;
01415         
01416         qmPCVecMap.insert(std::pair<Real,std::vector<int> >(
01417             currQMID, currPCSel ));
01418         
01419         numParsedPBDs++;
01420         delete customPCPDB;
01421     }
01422     
01423     delete pdbP;
01424     
01425     if (numParsedPBDs != qmNumGrps && simParams->qmCustomPCSel) {
01426         iout << iWARN << "The number of files provided for custom point "
01427         "charges is not equal to the number of QM groups!\n" << endi;
01428     }
01429     
01430     // Initializes an array with the number of Custom Point Charges per 
01431     // QM group.
01432     qmCustPCSizes = new int[qmNumGrps];
01433     for (int i=0; i<qmNumGrps; i++)
01434         qmCustPCSizes[i] = 0;
01435     
01436     qmTotCustPCs = 0;
01437     
01438     // Stores the size of each Custom PC vector in the array.
01439     // We may not have PCs for all QM groups.
01440     for (auto it = qmPCVecMap.begin(); it != qmPCVecMap.end(); it++) {
01441         qmTotCustPCs += (*it).second.size();
01442         int qmIndx = qmGrpIDMap[(*it).first];
01443         qmCustPCSizes[qmIndx] = (*it).second.size();
01444     }
01445     
01446     qmCustomPCIdxs = new int[qmTotCustPCs];
01447     
01448     if (simParams->qmCustomPCSel) {
01449         
01450         int auxIter = 0;
01451         for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
01452             
01453             Real currQMID = qmGrpID[grpIndx];
01454             
01455             iout << iINFO << "Loading " << qmPCVecMap[currQMID].size()
01456             << " custom point charges to QM Group " << grpIndx
01457             << " (ID: " << currQMID << ")\n" << endi;
01458             
01459             for (int pcIndx=0; pcIndx<qmPCVecMap[currQMID].size(); pcIndx++) {
01460                 qmCustomPCIdxs[auxIter] = qmPCVecMap[currQMID][pcIndx] ;
01461                 auxIter++;
01462             }
01463         }
01464     } 
01465     
01466     
01467     if (simParams->qmCSMD) {
01468         qmcSMD = simParams->qmCSMD;
01469         read_qm_csdm_file(qmGrpIDMap);
01470     }
01471     
01474     
01475     if (simParams->qmElecEmbed) {
01476         // Modifies Atom charges for Electrostatic Embedding.
01477         // QM atoms cannot have charges in the standard location, to keep
01478         // NAMD from calculating electrostatic interactions between QM and MM atoms.
01479         // We handle electrostatics ourselves in ComputeQM.C and in special
01480         // modifications for PME.
01481         for (int i=0; i<qmNumQMAtoms; i++) {
01482             qmAtmChrg[i] = atoms[qmAtmIndx[i]].charge;
01483             atoms[qmAtmIndx[i]].charge = 0;
01484         }
01485     }
01486     
01487     
01488     if ( simParams->extraBondsOn) {
01489         // Lifted from Molecule::build_extra_bonds
01490         
01491         StringList *file = cfgList->find("extraBondsFile");
01492         
01493         char err_msg[512];
01494         int a1,a2; float k, ref;
01495         
01496         for ( ; file; file = file->next ) {  // loop over files
01497             FILE *f = fopen(file->data,"r");
01498 //             if ( ! f ) {
01499 //               sprintf(err_msg, "UNABLE TO OPEN EXTRA BONDS FILE %s", file->data);
01500 //               NAMD_err(err_msg);
01501 //             } else {
01502 //               iout << iINFO << "READING EXTRA BONDS FILE " << file->data <<"\n"<<endi;
01503 //             }
01504             
01505             while ( 1 ) {
01506               char buffer[512];
01507               int ret_code;
01508               do {
01509                 ret_code = NAMD_read_line(f, buffer);
01510               } while ( (ret_code==0) && (NAMD_blank_string(buffer)) );
01511               if (ret_code!=0) break;
01512 
01513               char type[512];
01514               sscanf(buffer,"%s",type);
01515               
01516               int badline = 0;
01517               if ( ! strncasecmp(type,"bond",4) ) {
01518                 if ( sscanf(buffer, "%s %d %d %f %f %s",
01519                     type, &a1, &a2, &k, &ref, err_msg) != 5 ) badline = 1;
01520                 
01521                 // If an extra bond is defined between QM atoms, we make
01522                 // note so that it wont be deleted when we delete bonded 
01523                 // interactions between QM atoms.
01524                 if( qmAtomGroup[a1] > 0 && qmAtomGroup[a2]) {
01525                     Bond tmp;
01526                     tmp.bond_type = 0;
01527                     tmp.atom1 = a1;  tmp.atom2 = a2;
01528                     qmExtraBonds.add(tmp);
01529                 }
01530                 
01531                 
01532               }else if ( ! strncasecmp(type,"#",1) ) {
01533                 continue;  // comment
01534               }
01535               
01536             }
01537             fclose(f);
01538         }
01539     }
01540     
01541     return;
01542     
01543 #endif // MEM_OPT_VERSION
01544 }
01545 
01546 // Adapted from Molecule::delete_alch_bonded
01547 void Molecule::delete_qm_bonded(void)  {
01548 
01549 #ifdef MEM_OPT_VERSION
01550     NAMD_die("QMMM interface is not supported in memory optimized builds.");
01551 #else
01552     
01553     DebugM(3,"Cleaning QM bonds, angles, dihedrals, impropers and crossterms.\n")
01554     
01555     // Bonds
01556     Bond *nonQMBonds;
01557     nonQMBonds = new Bond[numBonds] ;
01558     int nonQMBondCount = 0;
01559     qmDroppedBonds = 0;
01560     for (int i = 0; i < numBonds; i++) {
01561         
01562         int part1 = qmAtomGroup[bonds[i].atom1];
01563         int part2 = qmAtomGroup[bonds[i].atom2];
01564         if (part1 > 0 && part2 > 0 ) {
01565             
01566             // If the user defined extra bonds, we check if the QM-QM bond is an "extra"
01567             // bond, and do not delete it.
01568             if (simParams->extraBondsOn) {
01569 //                 std::cout << "Checking Bond: " << bonds[i].atom1 << "," << bonds[i].atom2 << std::endl ;
01570                 
01571                 for (int ebi=0; ebi < qmExtraBonds.size() ; ebi++) {
01572                     
01573                     if( (qmExtraBonds[ebi].atom1 == bonds[i].atom1 || 
01574                         qmExtraBonds[ebi].atom1 == bonds[i].atom2) && 
01575                     (qmExtraBonds[ebi].atom2 == bonds[i].atom1 || 
01576                         qmExtraBonds[ebi].atom2 == bonds[i].atom2) ) {
01577 //                         std::cout << "This is an extra bond! We will keep it." << std::endl;
01578                         nonQMBonds[nonQMBondCount++] = bonds[i];
01579                         break;
01580                     }
01581                 }
01582             }
01583             
01584             qmDroppedBonds++;
01585         } else {
01586             // Just a simple sanity check.
01587             // If there are no QM bonds defined by the user but we find a
01588             // bond between a QM and an MM atom, error out.
01589             if ((part1 > 0 || part2 > 0) && qmNumBonds == 0) {
01590                 iout << iERROR << " Atoms " << bonds[i].atom1
01591                 << " and " << bonds[i].atom2 << " form a QM-MM bond!\n" << endi ;
01592                 NAMD_die("No QM-MM bonds were defined by the user, but one was found.");
01593             }
01594             nonQMBonds[nonQMBondCount++] = bonds[i];
01595         }
01596     }
01597     numBonds = nonQMBondCount;
01598     delete [] bonds;
01599     bonds = new Bond[numBonds];
01600     for (int i = 0; i < nonQMBondCount; i++) {
01601         bonds[i]=nonQMBonds[i];
01602     }
01603     delete [] nonQMBonds;
01604   numRealBonds = numBonds;
01605     
01606   // Angles
01607   Angle *nonQMAngles;
01608   nonQMAngles = new Angle[numAngles];
01609   int nonQMAngleCount = 0;
01610   qmDroppedAngles = 0;
01611   for (int i = 0; i < numAngles; i++) {
01612     int part1 = qmAtomGroup[angles[i].atom1];
01613     int part2 = qmAtomGroup[angles[i].atom2];
01614     int part3 = qmAtomGroup[angles[i].atom3];
01615     if (part1 > 0 && part2 > 0 && part3 > 0) {
01616       //CkPrintf("-----ANGLE ATOMS %i %i %i partitions %i %i %i\n",angles[i].atom1, angles[i].atom2, angles[i].atom3, part1, part2, part3);
01617       qmDroppedAngles++;
01618     }
01619     else {
01620       nonQMAngles[nonQMAngleCount++] = angles[i];
01621     }
01622   }
01623   numAngles = nonQMAngleCount;
01624   delete [] angles;
01625   angles = new Angle[numAngles];
01626   for (int i = 0; i < nonQMAngleCount; i++) {
01627     angles[i]=nonQMAngles[i];
01628   }
01629   delete [] nonQMAngles;
01630 
01631 
01632   // Dihedrals
01633   Dihedral *nonQMDihedrals;
01634   nonQMDihedrals = new Dihedral[numDihedrals];
01635   int nonQMDihedralCount = 0;
01636   qmDroppedDihedrals = 0;
01637   for (int i = 0; i < numDihedrals; i++) {
01638     int part1 = qmAtomGroup[dihedrals[i].atom1];
01639     int part2 = qmAtomGroup[dihedrals[i].atom2];
01640     int part3 = qmAtomGroup[dihedrals[i].atom3];
01641     int part4 = qmAtomGroup[dihedrals[i].atom4];
01642     if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
01643       //CkPrintf("-----i %i DIHEDRAL ATOMS %i %i %i %i partitions %i %i %i %i\n",i,dihedrals[i].atom1, dihedrals[i].atom2, dihedrals[i].atom3, dihedrals[i].atom4, part1, part2, part3,part4);
01644       qmDroppedDihedrals++;
01645     }
01646     else {
01647       nonQMDihedrals[nonQMDihedralCount++] = dihedrals[i];
01648     }
01649   }
01650   numDihedrals = nonQMDihedralCount;
01651   delete [] dihedrals;
01652   dihedrals = new Dihedral[numDihedrals];
01653   for (int i = 0; i < numDihedrals; i++) {
01654     dihedrals[i]=nonQMDihedrals[i];
01655   }
01656   delete [] nonQMDihedrals;
01657 
01658   // Impropers
01659   Improper *nonQMImpropers;
01660   nonQMImpropers = new Improper[numImpropers];
01661   int nonQMImproperCount = 0;
01662   qmDroppedImpropers = 0;
01663   for (int i = 0; i < numImpropers; i++) {
01664     int part1 = qmAtomGroup[impropers[i].atom1];
01665     int part2 = qmAtomGroup[impropers[i].atom2];
01666     int part3 = qmAtomGroup[impropers[i].atom3];
01667     int part4 = qmAtomGroup[impropers[i].atom4];
01668     if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0) {
01669       //CkPrintf("-----i %i IMPROPER ATOMS %i %i %i %i partitions %i %i %i %i\n",i,impropers[i].atom1, impropers[i].atom2, impropers[i].atom3, impropers[i].atom4, part1, part2, part3,part4);
01670       qmDroppedImpropers++;
01671     }
01672     else {
01673       nonQMImpropers[nonQMImproperCount++] = impropers[i];
01674     }
01675   }
01676   numImpropers = nonQMImproperCount;
01677   delete [] impropers;
01678   impropers = new Improper[numImpropers];
01679   for (int i = 0; i < numImpropers; i++) {
01680     impropers[i]=nonQMImpropers[i];
01681   }
01682   delete [] nonQMImpropers;
01683   
01684   // Crossterms 
01685   Crossterm *nonQMCrossterms;
01686   nonQMCrossterms = new Crossterm[numCrossterms];
01687   int nonQMCrosstermCount = 0;
01688   qmDroppedCrossterms = 0 ;
01689   for (int i = 0; i < numCrossterms; i++) {
01690     int part1 = qmAtomGroup[crossterms[i].atom1];
01691     int part2 = qmAtomGroup[crossterms[i].atom2];
01692     int part3 = qmAtomGroup[crossterms[i].atom3];
01693     int part4 = qmAtomGroup[crossterms[i].atom4];
01694     int part5 = qmAtomGroup[crossterms[i].atom5];
01695     int part6 = qmAtomGroup[crossterms[i].atom6];
01696     int part7 = qmAtomGroup[crossterms[i].atom7];
01697     int part8 = qmAtomGroup[crossterms[i].atom8];
01698     if (part1 > 0 && part2 > 0 && part3 > 0 && part4 > 0 &&
01699         part5 > 0 && part6 > 0 && part7 > 0 && part8 > 0) {
01700       //CkPrintf("-----i %i IMPROPER ATOMS %i %i %i %i partitions %i %i %i %i\n",i,impropers[i].atom1, impropers[i].atom2, impropers[i].atom3, impropers[i].atom4, part1, part2, part3,part4);
01701       qmDroppedCrossterms++;
01702     }
01703     else {
01704       nonQMCrossterms[nonQMCrosstermCount++] = crossterms[i];
01705     }
01706   }
01707   numCrossterms = nonQMCrosstermCount;
01708   delete [] crossterms;
01709   crossterms = new Crossterm[numCrossterms] ;
01710   for (int i=0; i < numCrossterms; i++){
01711       crossterms[i] = nonQMCrossterms[i] ;
01712   }
01713   delete [] nonQMCrossterms ;
01714   
01715   if (!CkMyPe()) {
01716       iout << iINFO << "The QM region will remove " << qmDroppedBonds << " bonds, " << 
01717           qmDroppedAngles << " angles, " << qmDroppedDihedrals << " dihedrals, "
01718           << qmDroppedImpropers << " impropers and " << qmDroppedCrossterms 
01719           << " crossterms.\n" << endi ;
01720   }
01721   
01722 #endif
01723 }
01724 
01725 typedef std::pair<int,int> cSMDPair ;
01726 
01727 void Molecule::read_qm_csdm_file(std::map<Real,int> &qmGrpIDMap)  {
01728     
01729     std::ifstream cSMDInputFile ;
01730     std::string Line("") ;
01731     
01732     iout << iINFO << "Reading QM cSMD configuration file: " << 
01733         simParams->qmCSMDFile << "\n" << endi ;
01734     
01735     cSMDInputFile.open( simParams->qmCSMDFile ) ;
01736     
01737     if (! cSMDInputFile.good())
01738     {
01739         NAMD_die("Configuration file for QM cSMD could not be read!") ;
01740     }
01741     
01742     
01743     // Index of Conditional SMD guides per group
01744     ResizeArray<ResizeArray<int> > cSMDindexV;
01745     cSMDindexV.resize(qmNumGrps);
01746     // Atom indices for Origin and Target of cSMD
01747     ResizeArray<cSMDPair> cSMDpairsV;
01748     // Spring constants for cSMD
01749     ResizeArray<Real> cSMDKsV;
01750     // Speed of movement of guide particles for cSMD.
01751     ResizeArray<Real> cSMDVelsV;
01752     // Distance cutoff for guide particles for cSMD.
01753     ResizeArray<Real> cSMDcoffsV;
01754     
01755     while( ! cSMDInputFile.eof() )
01756     {
01757         getline(cSMDInputFile, Line) ;
01758         
01759         if ( ! Line.length() )
01760             continue;
01761         
01762         if (Line.substr(0,1) == std::string("#") )
01763             continue;
01764         
01765         auto strVec = split( Line, " ");
01766         
01767         if (strVec.size() != 5 ) {
01768             iout << iERROR << "Format error in QM cSDM configuration file: " 
01769             << Line << "\n" << endi;
01770             NAMD_die("Error processing QM information.");
01771         }
01772         
01773         std::stringstream storConv ;
01774         Real convData ;
01775         
01776         cSMDpairsV.add( cSMDPair(0,0) );
01777         cSMDKsV.add(0);
01778         cSMDVelsV.add(0);
01779         cSMDcoffsV.add(0);
01780         
01781         for (int i=0; i < 5; i++ ) {
01782             storConv.clear() ;
01783             storConv << strVec[i] ;
01784             storConv >> convData;
01785             if (storConv.fail()) {
01786                 iout << iERROR << "Error parsing QM cSMD configuration file: " 
01787                 << convData << "\n" << endi;
01788                 NAMD_die("Error processing QM information.");
01789             }
01790             
01791             switch (i) {
01792             case 0:
01793                 cSMDpairsV[cSMDnumInst].first = convData ;
01794                 break;
01795             case 1:
01796                 cSMDpairsV[cSMDnumInst].second = convData ;
01797                 break;
01798             case 2:
01799                 cSMDKsV[cSMDnumInst] = convData ;
01800                 break;
01801             case 3:
01802                 cSMDVelsV[cSMDnumInst] = convData ;
01803                 break;
01804             case 4:
01805                 cSMDcoffsV[cSMDnumInst] = convData ;
01806                 break;
01807             }
01808         }
01809         
01810         // Sanity check
01811         if (cSMDpairsV[cSMDnumInst].first == cSMDpairsV[cSMDnumInst].second) {
01812             iout << iERROR << "Conditional SMD atoms must be different! We got " << 
01813             cSMDpairsV[cSMDnumInst].first << " and " << cSMDpairsV[cSMDnumInst].second << "\n" << endi;
01814             NAMD_die("Error processing QM information.") ;
01815         }
01816         
01817         // Sanity check
01818         if (qmAtomGroup[cSMDpairsV[cSMDnumInst].first] == 0 ||
01819             qmAtomGroup[cSMDpairsV[cSMDnumInst].second] == 0) {
01820             iout << iERROR << "Atoms " << cSMDpairsV[cSMDnumInst].first << " and " << 
01821             cSMDpairsV[cSMDnumInst].second << " MUST be assigned as QM atoms.\n" << endi;
01822             NAMD_die("Error processing QM information.") ;
01823         }
01824         
01825         // Sanity check
01826         if (qmAtomGroup[cSMDpairsV[cSMDnumInst].first] != 
01827             qmAtomGroup[cSMDpairsV[cSMDnumInst].second] ) {
01828             iout << iERROR << "Atoms in cSMD MUST be assigned to the same QM group.\n" << endi;
01829             NAMD_die("Error processing QM information.") ;
01830         }
01831         
01832         int grpIndx = qmGrpIDMap[qmAtomGroup[cSMDpairsV[cSMDnumInst].first]] ;
01833         
01834         iout << iINFO << "Adding cSMD data: (" << cSMDnumInst << ") " << 
01835         cSMDpairsV[cSMDnumInst].first << "," << cSMDpairsV[cSMDnumInst].second << "," << 
01836         cSMDKsV[cSMDnumInst] << "," << cSMDVelsV[cSMDnumInst] << "," << cSMDcoffsV[cSMDnumInst] <<
01837         " to QM Group " << grpIndx << "\n" << endi ;
01838         
01839         cSMDindexV[grpIndx].add(cSMDnumInst);
01840         
01841         cSMDnumInst++;
01842     }
01843     
01844     cSMDindex = new int*[qmNumGrps];
01845     cSMDindxLen = new int[qmNumGrps];
01846     
01847     for (int i=0; i<qmNumGrps; i++) {
01848             cSMDindex[i] = new int[cSMDindexV[i].size()];
01849             cSMDindxLen[i] = cSMDindexV[i].size();
01850     
01851             for (int j=0; j<cSMDindxLen[i]; j++)
01852                     cSMDindex[i][j] = cSMDindexV[i][j] ;
01853     }
01854     
01855     cSMDpairs = new int*[cSMDnumInst];
01856     for (int i=0; i<cSMDnumInst; i++)
01857             cSMDpairs[i] = new int[2];
01858     cSMDKs  = new Real[cSMDnumInst]; 
01859     cSMDVels = new Real[cSMDnumInst];
01860     cSMDcoffs = new Real[cSMDnumInst];
01861 
01862     for (int i=0; i<cSMDnumInst; i++) {
01863             cSMDpairs[i][0] = cSMDpairsV[i].first;
01864             cSMDpairs[i][1] = cSMDpairsV[i].second;
01865     
01866             cSMDKs[i] = cSMDKsV[i];
01867             cSMDVels[i] = cSMDVelsV[i];
01868             cSMDcoffs[i] = cSMDcoffsV[i];
01869     }
01870     
01871 }
01872 
01873 
01874 
01875 
01876 

Generated on Fri Sep 22 01:17:13 2017 for NAMD by  doxygen 1.4.7