Fix QM/MM memory allocation error
[namd.git] / src / MoleculeQM.C
1 /**
2 ***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
3 ***  The Board of Trustees of the University of Illinois.
4 ***  All rights reserved.
5 **/
6
7 #ifdef DEBUG_QM
8   #define DEBUGM
9 #endif
10
11 #include <stdio.h>
12 #include <string.h>
13 #include <stdlib.h>
14 #include <ctype.h>
15 #include "ResizeArray.h"
16 #include "InfoStream.h"
17 #include "Molecule.h"
18 #include "strlib.h"
19 #include "MStream.h"
20 #include "Communicate.h"
21 #include "Node.h"
22 #include "ObjectArena.h"
23 #include "Parameters.h"
24 #include "PDB.h"
25 #include "SimParameters.h"
26 #include "Hydrogen.h"
27 #include "UniqueSetIter.h"
28 #include "charm++.h"
29
30 #define MIN_DEBUG_LEVEL 3
31 //#define DEBUGM
32 #include "Debug.h"
33 #include "CompressPsf.h"
34 #include "ParallelIOMgr.h"
35 #include <deque>
36 #include <algorithm>
37
38 #ifndef CODE_REDUNDANT
39 #define CODE_REDUNDANT 0
40 #endif
41
42 #include <string>
43 #include <sstream>
44 #include <fstream>
45
46 class qmSolvData {
47 public:
48     char segName[5];
49     int resID, begAtmID, size;
50     std::vector<int> atmIDs ;
51     Real qmGrpID ;
52     
53     qmSolvData() : resID(-1), begAtmID(-1), size(0) {}
54     qmSolvData(int newResID, int newBegAtm, int newSize) {
55         resID = newResID;
56         begAtmID = newBegAtm;
57         size = newSize;
58     }
59     qmSolvData(int newResID, int newBegAtm, int newSize, 
60                const char* newSegName, Real newQMID) {
61         resID = newResID;
62         begAtmID = newBegAtm;
63         size = newSize;
64         strncpy(segName, newSegName,5);
65         atmIDs.push_back(newBegAtm) ;
66         qmGrpID = newQMID;
67     }
68     
69     bool operator <(const qmSolvData& ref) {return begAtmID < ref.begAtmID;}
70     bool operator ==(const qmSolvData& ref) {return begAtmID == ref.begAtmID;}
71 };
72
73 std::vector<std::string> split(const std::string &text, std::string delimiter) {
74     
75     std::vector<std::string> tokens;
76     std::size_t start = 0, end = 0;
77     
78     while ((end = text.find(delimiter, start)) != std::string::npos) {
79         
80         std::string temp = text.substr(start, end - start);
81         
82         if (! temp.empty()) tokens.push_back(temp);
83         
84         start = end + delimiter.length();
85     }
86     
87     // Gets last item
88     std::string temp = text.substr(start);
89     if (! temp.empty()) tokens.push_back(temp);
90     
91     return tokens;
92 }
93
94 struct refSelStr {
95     std::string segid;
96     int resid;
97     
98     refSelStr() {} ;
99     refSelStr(std::string newSegID, int newResid) {
100         segid = newSegID;
101         resid = newResid;
102     }
103 } ;
104
105 typedef std::vector<refSelStr> refSelStrVec ;
106 typedef std::map<Real, refSelStrVec> refSelStrMap ;
107 typedef std::pair<Real,refSelStrVec> refSelStrPair ;
108
109 void Molecule::prepare_qm(const char *pdbFileName, 
110                                Parameters *params, ConfigList *cfgList) {
111 #ifdef MEM_OPT_VERSION
112     NAMD_die("QMMM interface is not supported in memory optimized builds.");
113 #else
114     
115     PDB *pdbP;      //  Pointer to PDB object to use
116     
117     qmNumQMAtoms = 0;
118     
119     qmNumGrps = 0 ;
120     
121     iout << iINFO << "Using the following PDB file for QM parameters: " << 
122     pdbFileName << "\n" << endi;
123     if (pdbFileName)
124         pdbP = new PDB(pdbFileName);
125     else
126         iout << "pdbFile name not defined!\n\n" << endi ;
127     
128     if (pdbP->num_atoms() != numAtoms)
129     {
130        NAMD_die("Number of atoms in QM parameter PDB file doesn't match coordinate PDB");
131     }
132     
133     qmElementArray = new char*[numAtoms] ;
134     
135     qmAtomGroup = new Real[numAtoms] ;
136     
137     BigReal bondVal;
138     Real atmGrp;
139     
140     std::set<Real> atmGrpSet;
141     
142     std::vector<Real> grpChrgVec;
143     
144     // Value in the bond column fro QM-MM bonds.
145     std::vector<BigReal> qmBondValVec;
146     // Identifiers of different QM regions
147     std::vector<Real> qmGroupIDsVec;
148     // This maps the qm group ID with the serail index of this group in 
149     // various arrays.
150     std::map<Real,int> qmGrpIDMap ;
151     
152     std::vector<int> qmAtmIndxVec, qmGrpSizeVec;
153     
154     // Used to parse user selection in different places.
155     std::vector<std::string> strVec;
156     
157     qmNoPC = simParams->qmNoPC;
158     
159     // We set to zero by default, So there is no need for extra processing.
160     qmPCFreq = 0;
161     if (simParams->qmPCSelFreq > 1)
162         qmPCFreq = simParams->qmPCSelFreq;
163     
164     
165     ///////////////////////////////
166     /// Prepares Live Solvent Selection
167     
168     
169     qmLSSTotalNumAtms = 0;
170     SortedArray< qmSolvData> lssGrpRes;
171     std::map<Real,std::vector<int> > lssGrpRefIDs;
172     refSelStrMap lssRefUsrSel;
173     int totalNumRefAtms = 0;
174     int lssClassicResIndx = 0 ;
175     int lssCurrClassResID = -1 ;
176     char lssCurrClassResSegID[5];
177     if (simParams->qmLSSOn) {
178         DebugM(4, "LSS is ON! Processing QM solvent.\n") ;
179         
180         if (resLookup == NULL)
181             NAMD_die("We need residue data to conduct LSS.");
182          
183         if (simParams->qmLSSMode == QMLSSMODECOM) {
184             
185             StringList *current = cfgList->find("QMLSSRef");
186             for ( ; current; current = current->next ) {
187                 
188                 strVec = split( std::string(current->data) , " ");
189                 
190                 if (strVec.size() != 3 ) {
191                     iout << iERROR << "Format error in QM LSS size: " 
192                     << current->data
193                     << "\n" << endi;
194                     NAMD_die("Error processing QM information.");
195                 }
196                 
197                 std::stringstream storConv ;
198                 
199                 storConv << strVec[0] ;
200                 Real grpID ;
201                 storConv >> grpID;
202                 if (storConv.fail()) {
203                     iout << iERROR << "Error parsing QMLSSRef selection: " 
204                     << current->data
205                     << "\n" << endi;
206                     NAMD_die("Error processing QM information.");
207                 }
208                 
209                 std::string segName = strVec[1].substr(0,4);
210                 
211                 storConv.clear() ;
212                 storConv << strVec[2];
213                 int resID ;
214                 storConv >> resID;
215                 if (storConv.fail()) {
216                     iout << iERROR << "Error parsing QMLSSRef selection: " 
217                     << current->data
218                     << "\n" << endi;
219                     NAMD_die("Error processing QM information.");
220                 }
221                 
222                 auto it = lssRefUsrSel.find(grpID) ;
223                 if (it == lssRefUsrSel.end())
224                     lssRefUsrSel.insert(refSelStrPair(grpID,refSelStrVec()));
225                 
226                 lssRefUsrSel[grpID].push_back(refSelStr(segName, resID));
227             }
228             
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 ;
235                 }
236             }
237         }
238     }
239     
240     
241     
242     ///////////////////////////////
243     /// Data gathering from PDB to find QM atom and bond info
244     
245     
246     for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
247         
248         // If we are looking for QM-MM bonds, then we need to store extra info.
249         if (simParams->qmBondOn) {
250             
251             // We store both the qm group and the bond value
252             if ( strcmp(simParams->qmColumn,"beta") == 0 ){
253                 atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
254                 
255                 bondVal = pdbP->atom(atmInd)->occupancy() ;
256             }
257             else {
258                 atmGrp = pdbP->atom(atmInd)->occupancy() ;
259                 
260                 bondVal = pdbP->atom(atmInd)->temperaturefactor() ;
261             }
262             
263             qmBondValVec.push_back(bondVal);
264         }
265         else {
266             
267             if ( strcmp(simParams->qmColumn,"beta") == 0 ){
268                 atmGrp = pdbP->atom(atmInd)->temperaturefactor() ;
269             }
270             else {
271                 atmGrp = pdbP->atom(atmInd)->occupancy() ;
272             }
273         }
274         
275         // We store all atom QM Group IDs in an array for 
276         // future transport to all PEs.
277         qmAtomGroup[atmInd] = atmGrp;
278         
279         // We store all atom's elements for quick access in the QM code.
280         // Elements are used to tell the QM code the atom type.
281         qmElementArray[atmInd] = new char[3];
282         strncpy(qmElementArray[atmInd],pdbP->atom(atmInd)->element(),3);
283         
284         // For QM atoms, we keep global atom indices and parse different QM Group
285         // IDs, keeping track of each groups size
286         if (atmGrp > 0){
287             
288             if (simParams->fixedAtomsOn) {
289                 if ( fixedAtomFlags[atmInd] == 1 ) {
290                     iout << iERROR << "QM atom cannot be fixed in space!\n" 
291                     << endi;
292                     NAMD_die("Error processing QM information.");
293                 }
294             }
295             
296             // Changes the VdW type of QM atoms.
297             // This is may be used to counter balance the overpolarization 
298             // that QM atoms sufer.
299             if (simParams->qmVDW) {
300                 // Creating a new type string
301                 std::string qmType(qmElementArray[atmInd]) ;
302                 // Erases empty spaces
303                 qmType.erase(
304                         std::remove_if(qmType.begin(), qmType.end(), isspace ),
305                     qmType.end());
306                 // pre-pends a "q" to the element name.
307                 qmType = std::string("q") + qmType;
308             
309 //                 iout << "QM VdW type: " << atoms[atmInd].vdw_type 
310 //                 << " atom type: " << atomNames[atmInd].atomtype 
311 //                 << " new type "<< qmType << "\n" << endi;
312                 
313                 /*  Look up the vdw constants for this atom    */
314                 // I am casting a non-const here because the function doesn't actually
315                 // change the string value, but it doesn't ask for a const char* as 
316                 // an argument.
317                 params->assign_vdw_index(const_cast<char*>( qmType.c_str()), 
318                                          &(atoms[atmInd]));
319                 
320 //                 iout << "----> new VdW type: " << atoms[atmInd].vdw_type << "\n" << endi;
321             }
322             
323             // Indexes all global indices of QM atoms.
324             qmAtmIndxVec.push_back(atmInd);
325             
326             auto retVal = atmGrpSet.insert(atmGrp);
327             
328             // If a new QM group is found, store the reverse reference in a map
329             // and increase the total count.
330             if (retVal.second) {
331                 // This mak makes the reverse identification from the group ID
332                 // to the sequential order in which each group was first found.
333                 qmGrpIDMap.insert(std::pair<BigReal,int>(atmGrp, qmNumGrps));
334                 
335                 // This vector keeps track of the group ID for each group
336                 qmGroupIDsVec.push_back(atmGrp);
337                 
338                 // This counter is used to keep track of the sequential order in
339                 // which QM groups were first seen in the reference PDB file.
340                 qmNumGrps++ ;
341                 
342                 // If this is a new QM group, initialize a new variable in the 
343                 // vector to keep track of the number of atoms in this group.
344                 qmGrpSizeVec.push_back(1);
345                 
346                 // For each new QM group, a new entry in the total charge
347                 // vector is created
348                 grpChrgVec.push_back( atoms[atmInd].charge );
349             }
350             else {
351                 // If the QM group has already been seen, increase the count
352                 // of atoms in that group.
353                 qmGrpSizeVec[ qmGrpIDMap[atmGrp] ]++ ;
354                 
355                 // If the QM group exists, adds current atom charge to total charge.
356                 grpChrgVec[ qmGrpIDMap[atmGrp] ] += atoms[atmInd].charge;
357             }
358             
359             // In case we are using live solvent selection
360             if(simParams->qmLSSOn) {
361                 
362                 int resID = pdbP->atom(atmInd)->residueseq();
363                 char segName[5];
364                 strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
365                 
366                 int resSize = get_residue_size(segName,resID);
367                 
368                 int i =-1, end =-1;
369                 
370                 resLookup->lookup(segName,resID,&i, &end);
371                 
372                 if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
373                            simParams->qmLSSResname, 4) == 0) {
374                     
375                     // We try to insert every residue from every atom
376                     qmSolvData *resP = lssGrpRes.find(qmSolvData(resID, i, resSize));
377                     
378                     if (resP != NULL) {
379                         resP->atmIDs.push_back(atmInd) ;
380 //                         DebugM(3, "Existing residue " << resP->resID 
381 //                         << " from segID " << resP->segName
382 //                         << " received atom "
383 //                         << atmInd << "\n" );
384                     }
385                     else {
386                         lssGrpRes.insert(qmSolvData(resID, i, resSize, segName, atmGrp));
387 //                         DebugM(3, lssGrpRes.size() << ") Inserted new residue: " 
388 //                         << resID << " from segID " << segName
389 //                         << " with atom "
390 //                         << i << "\n" ) ;
391                     }
392                 }
393                 else {
394                     // We store the QM atoms, per group, which are NOT part of
395                     // solvent molecules.
396                     
397                     // Checks if we have a vector for this QM group.
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>() ) );
402                     }
403                     
404                     switch (simParams->qmLSSMode)
405                     {
406                     
407                     case QMLSSMODECOM:
408                     {
409                         // If we are using COM selection, checks if the atom
410                         // is part of the user selected residues
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(),
414                                         segName,5) == 0 ) {
415                                 
416                                 lssGrpRefIDs[atmGrp].push_back(atmInd);
417                                 totalNumRefAtms++;
418                             }
419                         }
420                         
421                     } break;
422                     
423                     case QMLSSMODEDIST:
424                     {
425                     
426                         lssGrpRefIDs[atmGrp].push_back(atmInd);
427                         totalNumRefAtms++;
428                     
429                     } break;
430                     
431                     }
432                 }
433                 
434             }
435             
436         }
437         else if (atmGrp == 0) {
438             
439             if(simParams->qmLSSOn) {
440                 
441                 if ( strncasecmp(pdbP->atom(atmInd)->residuename(),
442                            simParams->qmLSSResname, 4) == 0) {
443                     
444                     int resID = pdbP->atom(atmInd)->residueseq();
445                     char segName[5];
446                     strncpy(segName, pdbP->atom(atmInd)->segmentname(),5);
447                     
448                     if (lssCurrClassResID < 0) {
449                         lssCurrClassResID = resID ;
450                         strncpy(lssCurrClassResSegID, segName,5);
451                         lssClassicResIndx = 0;
452                     }
453                     else if (lssCurrClassResID != resID ||
454                             strcmp(lssCurrClassResSegID, segName) != 0 ) {
455                         lssCurrClassResID = resID ;
456                         strncpy(lssCurrClassResSegID, segName,5);
457                         lssClassicResIndx++;
458                     }
459                     
460                     qmClassicSolv.insert(std::pair<int,int>(atmInd,lssClassicResIndx));
461                     
462                 }
463             }
464         }
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.");
468         }
469     }
470     
471     // Sanity check
472     if (simParams->qmLSSOn) {
473         if (lssGrpRes.size() == 0)
474             NAMD_die("LSS was activated but there are no QM solvent molecules!\n");
475         
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.");
482             }
483         }
484         
485         DebugM(3,"We found " << lssClassicResIndx << " classical solvent molecules totalling "
486         << qmClassicSolv.size() << " atoms.\n" );
487         
488     }
489     
490     qmNumQMAtoms = qmAtmIndxVec.size();
491     
492     if (qmNumQMAtoms == 0)
493         NAMD_die("No QM atoms were found in this QM simulation!") ;
494     
495     iout << iINFO << "Number of QM atoms (excluding Dummy atoms): " << 
496     qmNumQMAtoms << "\n" << endi;
497     
498     qmAtmChrg = new Real[qmNumQMAtoms] ;
499     qmAtmIndx = new int[qmNumQMAtoms] ;
500     for (int i=0; i<qmNumQMAtoms; i++) {
501         // qmAtmChrg gets initialized with the PSF charges at the end of this
502         // function, but values may change as QM steps are taken.
503         qmAtmChrg[i] = 0;  
504         qmAtmIndx[i] = qmAtmIndxVec[i] ;
505     }
506     
507     // This map relates the QM group index with a vector of pairs
508     // of bonded MM-QM atoms (with the bonded atoms ins this order: 
509     // MM first, QM second).
510     std::map<int, std::vector<std::pair<int,int> > > qmGrpIDBonds ;
511     int bondCounter = 0;
512     if (simParams->qmBondOn) {
513         
514         // Checks all bonds for QM-MM pairs.
515         // Makes sanity checks against QM-QM pairs and MM-MM pairs that
516         // are flagged by the user to be bonds between QM and MM regions.
517         // QM-QM bonds will be removed in another function.
518         for (int bndIt = 0; bndIt < numBonds; bndIt++) {
519             
520             bond curr = bonds[bndIt] ;
521             
522             // In case either atom has a non-zero
523             if ( qmBondValVec[curr.atom1] != 0 &&
524                  qmBondValVec[curr.atom2] != 0 ) {
525                 
526                 // Sanity checks (1 of 2)
527                 if (qmAtomGroup[curr.atom1] != 0 &&
528                     qmAtomGroup[curr.atom2] != 0) {
529                     iout << iERROR << "Atoms " << curr.atom1 << " and " << 
530                     curr.atom2 << " are assigned as QM atoms.\n" << endi;
531                     NAMD_die("Error in QM-MM bond assignment.") ;
532                 }
533                 
534                 // Sanity checks (2 of 2)
535                 if (qmAtomGroup[curr.atom1] == 0 &&
536                     qmAtomGroup[curr.atom2] == 0) {
537                     iout << iERROR << "Atoms " << curr.atom1 << " and " << 
538                     curr.atom2 << " are assigned as MM atoms.\n" << endi;
539                     NAMD_die("Error in QM-MM bond assignment.") ;
540                 }
541                 
542                 int currGrpID = 0;
543                 std::pair<int,int> newPair(0,0);
544                 
545                 // We create a QM-MM pair with the MM atom first
546                 if (qmAtomGroup[curr.atom1] != 0) {
547                     newPair = std::pair<int,int>(curr.atom2, curr.atom1) ;
548                     currGrpID = qmAtomGroup[curr.atom1];
549                 } else {
550                     newPair = std::pair<int,int>(curr.atom1, curr.atom2) ;
551                     currGrpID = qmAtomGroup[curr.atom2];
552                 } 
553                 
554                 int grpIndx = qmGrpIDMap[currGrpID] ;
555                 
556                 // Stores bonds in vectors key-ed by the QM group ID of the QM atom.
557                 auto retIter = qmGrpIDBonds.find(grpIndx) ;
558                 // In case thi QM-MM bonds belong to a QM group we have not seen 
559                 // yet, stores a new vector in the map.
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> >() ) );
563                 }
564                 
565                 qmGrpIDBonds[grpIndx].push_back( newPair );
566                 
567                 bondCounter++ ;
568             }
569             
570             
571         }
572         
573 //         iout << "Finished processing QM-MM bonds.\n" << endi ;
574         
575         if (bondCounter == 0)
576             iout << iWARN << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
577         else
578             iout << iINFO << "We found " << bondCounter << " QM-MM bonds.\n" << endi ;
579         
580     }
581     
582     // Initializes several arrays used to setup the QM simulation.
583     qmNumBonds = bondCounter ;
584     
585     qmMMBond = new int*[qmNumBonds];
586     
587     qmDummyBondVal = new BigReal[qmNumBonds];
588     
589     qmMMChargeTarget = new int*[qmNumBonds] ;
590     qmMMNumTargs = new int[qmNumBonds] ;
591     
592     qmDummyElement = new char*[qmNumBonds] ;
593     
594     
595     qmNumGrps = atmGrpSet.size();
596     
597     qmGrpSizes = new int[qmNumGrps];
598     
599     qmGrpID = new Real[qmNumGrps];
600     
601     qmGrpChrg = new Real[qmNumGrps];
602     
603     qmGrpMult = new Real[qmNumGrps] ;
604     
605     qmGrpNumBonds = new int[qmNumGrps];
606     
607     qmGrpBonds = new int*[qmNumGrps];
608     qmMMBondedIndx = new int*[qmNumGrps];
609     
610     
611     ///////////////////////////////
612     /// Multiplicity of each QM group
613     
614     
615     // We first initialize the multiplicity vector with 
616     // default values, then populate it with user defined values.
617     for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
618         qmGrpMult[grpIter] = 0;
619     }
620     
621     int multCount = 0;
622     StringList *current = cfgList->find("QMMult");
623     for ( ; current; current = current->next ) {
624         
625         auto strVec = split( std::string(current->data) , " ");
626         
627         if (strVec.size() != 2 ) {
628             iout << iERROR << "Format error in QM Multiplicity string: " 
629             << current->data
630             << "\n" << endi;
631             NAMD_die("Error processing QM information.");
632         }
633         
634         std::stringstream storConv ;
635         
636         storConv << strVec[0] ;
637         Real grpID ;
638         storConv >> grpID;
639         if (storConv.fail()) {
640             iout << iERROR << "Error parsing QMMult selection: " 
641             << current->data
642             << "\n" << endi;
643             NAMD_die("Error processing QM information.");
644         }
645         
646         storConv.clear() ;
647         storConv << strVec[1];
648         Real multiplicity ;
649         storConv >> multiplicity;
650         if (storConv.fail()) {
651             iout << iERROR << "Error parsing QMMult selection: " 
652             << current->data
653             << "\n" << endi;
654             NAMD_die("Error processing QM information.");
655         }
656         
657         auto it = qmGrpIDMap.find(grpID);
658         
659         if (it == qmGrpIDMap.end()) {
660             iout << iERROR << "Error parsing QMMult selection. Could not find QM group ID: " 
661             << grpID
662             << "\n" << endi;
663             NAMD_die("Error processing QM information.");
664         }
665         else {
666             iout << iINFO << "Applying user defined multiplicity "
667             << multiplicity << " to QM group ID " << grpID << "\n" << endi;
668             qmGrpMult[it->second] = multiplicity;
669         }
670         
671         multCount++;
672     }
673     
674     if (multCount != qmNumGrps && simParams->qmFormat == QMFormatORCA) {
675         NAMD_die("ORCA neds multiplicity values for all QM regions.");
676     }
677     
678     
679     ///////////////////////////////
680     /// Charge of each QM group
681     
682     
683     for (size_t grpIndx = 0; grpIndx < qmGrpSizeVec.size(); grpIndx++) {
684         
685         bool nonInteger = true;
686         if ((fabsf(roundf(grpChrgVec[grpIndx]) - grpChrgVec[grpIndx]) <= 0.001f) ) {
687             grpChrgVec[grpIndx] = roundf(grpChrgVec[grpIndx]) ;
688             nonInteger = false;
689         }
690         
691         iout << iINFO << grpIndx + 1 << ") Group ID: " << qmGroupIDsVec[grpIndx]
692         << " ; Group size: " << qmGrpSizeVec[grpIndx] << " atoms"
693         << " ; Total PSF charge: " << grpChrgVec[grpIndx] << "\n" << endi ;
694         
695         if (nonInteger && simParams->PMEOn)
696             NAMD_die("QM atoms do not add up to a whole charge, which is needed for PME.") ;
697     }
698     
699     int chrgCount = 0;
700     current = cfgList->find("QMCharge");
701     for ( ; current; current = current->next ) {
702         
703         auto strVec = split( std::string(current->data) , " ");
704         
705         if (strVec.size() != 2 ) {
706             iout << iERROR << "Format error in QM Charge string: " 
707             << current->data
708             << "\n" << endi;
709             NAMD_die("Error processing QM information.");
710         }
711         
712         std::stringstream storConv ;
713         
714         storConv << strVec[0] ;
715         Real grpID ;
716         storConv >> grpID;
717         if (storConv.fail()) {
718             iout << iERROR << "Error parsing QMCharge selection: " 
719             << current->data
720             << "\n" << endi;
721             NAMD_die("Error processing QM information.");
722         }
723         
724         storConv.clear() ;
725         storConv << strVec[1];
726         Real charge ;
727         storConv >> charge;
728         if (storConv.fail()) {
729             iout << iERROR << "Error parsing QMCharge selection: " 
730             << current->data
731             << "\n" << endi;
732             NAMD_die("Error processing QM information.");
733         }
734         
735         auto it = qmGrpIDMap.find(grpID);
736         
737         if (it == qmGrpIDMap.end()) {
738             iout << iERROR << "Error parsing QMCharge selection. Could not find QM group ID: " 
739             << grpID
740             << "\n" << endi;
741             NAMD_die("Error processing QM information.");
742         }
743         else {
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;
747         }
748         
749         chrgCount++;
750     }
751     
752     simParams->qmMOPACAddConfigChrg = false;
753     // Checks if QM group charges can be defined for MOPAC.
754     // Since charge can be defined in QMConfigLine, we need extra logic.
755     if (simParams->qmFormat == QMFormatMOPAC) {
756         
757         // Checks if group charge was supplied in config line for MOPAC.
758         std::string::size_type chrgLoc = std::string::npos ;
759         
760         // This will hold a sting with the first user supplied configuration line for MOPAC.
761         std::string configLine(cfgList->find("QMConfigLine")->data) ;
762         chrgLoc = configLine.find("CHARGE") ;
763         
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;
768         }
769         else if ( chrgLoc == std::string::npos && (simParams->qmChrgFromPSF || chrgCount == qmNumGrps) ) {
770         // If no charge was defined in the configuration line, gets from PSF and/or 
771         // from user defined charges (through the QMCharge keyword).
772             simParams->qmMOPACAddConfigChrg = true;
773         }
774         else
775         {
776         // If we could nont find a charge definition in the config line AND
777         // no specific charge was selected for each QM region through QMCharge AND
778         // the QMChargeFromPSF was not turned ON, the we stop NAMD and scream at the user.
779 //         if ( chrgLoc == std::string::npos && (chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) 
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.");
783         }
784         
785     }
786     
787     if (simParams->qmFormat == QMFormatORCA ) {
788         if ((chrgCount != qmNumGrps ) && !simParams->qmChrgFromPSF) {
789             // If we are not supposed to get charges from the PSF, and not 
790             // enough charges were set to cover all QM regions, 
791             // we stop NAMD and scream at the user.
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.");
795         }
796     }
797     
798     // If mechanichal embedding was requested but we have QM-MM bonds, we need 
799     // to send extra info to ComputeQM to preserve calculation speed.
800     if (qmNumBonds > 0 && qmNoPC) {
801         qmMeNumBonds = qmNumBonds;
802         qmMeMMindx = new int[qmMeNumBonds] ;
803         qmMeQMGrp = new Real[qmMeNumBonds] ;
804     }
805     else {
806         qmMeNumBonds = 0 ;
807         qmMeMMindx = NULL;
808         qmMeQMGrp = NULL;
809     }
810     
811     
812     ///////////////////////////////
813     /// Populate arrays that are used throughout the the calculations.
814     
815     
816     bondCounter = 0;
817     for (int grpIter = 0; grpIter < qmNumGrps; grpIter++) {
818         
819         qmGrpID[grpIter] = qmGroupIDsVec[grpIter] ;
820         
821         qmGrpSizes[grpIter] = qmGrpSizeVec[grpIter] ;
822         
823         qmGrpChrg[grpIter] = grpChrgVec[grpIter];
824         
825 //         iout << "Loaded " << qmGrpSizes[grpIter] << " atoms to this group.\n" << endi;
826         
827         int currNumbBonds = qmGrpIDBonds[grpIter].size() ;
828         
829         // Assigns the number of bonds that the current QM group has.
830         qmGrpNumBonds[grpIter] = currNumbBonds;
831         
832         if (currNumbBonds > 0) {
833             
834             qmGrpBonds[grpIter] = new int[currNumbBonds];
835             qmMMBondedIndx[grpIter] = new int[currNumbBonds];
836             
837             for (int bndIter = 0; bndIter<currNumbBonds; bndIter++) {
838                 
839                 // Adds the bonds to the overall sequential list.
840                 qmMMBond[bondCounter] = new int[2] ;
841                 qmMMBond[bondCounter][0] = qmGrpIDBonds[grpIter][bndIter].first ;
842                 qmMMBond[bondCounter][1] = qmGrpIDBonds[grpIter][bndIter].second ;
843                 
844                 // For the current QM group, and the current bond, gets the bond index.
845                 qmGrpBonds[grpIter][bndIter] =  bondCounter;
846                 
847                 // For the current QM group, and the current bond, gets the MM atom.
848                 qmMMBondedIndx[grpIter][bndIter] = qmGrpIDBonds[grpIter][bndIter].first ;
849                 
850                 // Assign the default value of dummy element
851                 qmDummyElement[bondCounter] = new char[3];
852                 strcpy(qmDummyElement[bondCounter],"H\0");
853                 
854                 // Determines the distance that will separate the new Dummy atom
855                 // and the Qm atom to which it will be bound.
856                 bondVal = 0;
857                 if (simParams->qmBondDist) {
858                     if (qmBondValVec[qmMMBond[bondCounter][0]] != 
859                         qmBondValVec[qmMMBond[bondCounter][1]]
860                     ) {
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]
863                         << ".\n" << endi ;
864                         NAMD_die("Error in QM-MM bond processing.");
865                     }
866                     
867                     bondVal = qmBondValVec[qmMMBond[bondCounter][0]] ;
868                 } 
869                 else {
870                     
871                     if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"C") == 0 ) {
872                         bondVal = 1.09 ;
873                     }
874                     else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"O") == 0 ) {
875                         bondVal = 0.98 ;
876                     }
877                     else if (strcmp(qmElementArray[qmMMBond[bondCounter][1]],"N") == 0 ) {
878                         bondVal = 0.99 ;
879                     }
880                     else {
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]] 
884                         <<  ".\n" << endi ;
885                         NAMD_die("Error in QM-MM bond processing.");
886                     }
887                     
888                 }
889                 
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] << ")"
894                 << "\n" << endi ;
895                 
896                 qmDummyBondVal[bondCounter] = bondVal;
897                 
898                 // In case we are preparing for a mechanical embedding simulation
899                 // with no point charges, populate the following vectors
900                 if (qmMeNumBonds > 0) {
901                     qmMeMMindx[bondCounter] = qmMMBond[bondCounter][0];
902                     qmMeQMGrp[bondCounter] = qmGrpID[grpIter];
903                 }
904                 
905                 bondCounter++ ;
906             }
907         }
908         
909     }
910     
911     
912     ///////////////////////////////
913     /// Overides Link Atom element with user selection.
914     
915     current = NULL;
916     if (qmNumBonds > 0)
917         current = cfgList->find("QMLinkElement");
918     
919     int numParsedLinkElem = 0;
920     for ( ; current != NULL; current = current->next ) {
921         
922         DebugM(3,"Parsing link atom element: " << current->data << "\n" );
923         
924         strVec = split( std::string(current->data) , " ");
925         
926         // We need the two atoms that compose the QM-MM bonds and 
927         // then the element.
928         if (strVec.size() != 3 && qmNumBonds > 1) {
929             iout << iERROR << "Format error in QM link atom element selection: " 
930             << current->data
931             << "\n" << endi;
932             NAMD_die("Error processing QM information.");
933         }
934         
935         // If there is only one QM-MM bond, we can accept the element only.
936         if (strVec.size() != 1 && strVec.size() != 3 && qmNumBonds == 1) {
937             iout << iERROR << "Format error in QM link atom element selection: " 
938             << current->data
939             << "\n" << endi;
940             NAMD_die("Error processing QM information.");
941         }
942         
943         std::stringstream storConv ;
944         int bondAtom1, bondAtom2;
945         std::string linkElement ;
946         
947         if (strVec.size() == 1) {
948             linkElement = strVec[0].substr(0,2);
949         }
950         else if (strVec.size() == 3) {
951             
952             storConv << strVec[0] ;
953             storConv >> bondAtom1;
954             if (storConv.fail()) {
955                 iout << iERROR << "Error parsing link atom element selection: " 
956                 << current->data
957                 << "\n" << endi;
958                 NAMD_die("Error processing QM information.");
959             }
960             
961             storConv.clear() ;
962             storConv << strVec[1];
963             storConv >> bondAtom2;
964             if (storConv.fail()) {
965                 iout << iERROR << "Error parsing link atom element selection: " 
966                 << current->data
967                 << "\n" << endi;
968                 NAMD_die("Error processing QM information.");
969             }
970             
971             linkElement = strVec[2].substr(0,2);
972         }
973         
974         numParsedLinkElem++;
975         
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.");
980         }
981         
982         int bondIter;
983         
984         if (strVec.size() == 1) {
985             bondIter = 0;
986         }
987         else if (strVec.size() == 3) {
988             
989             Bool foundBond = false;
990             
991             for (bondIter=0; bondIter<qmNumBonds; bondIter++) {
992                 
993                 if ( (qmMMBond[bondIter][0] == bondAtom1 &&
994                     qmMMBond[bondIter][1] == bondAtom2 ) ||
995                     (qmMMBond[bondIter][0] == bondAtom2 &&
996                     qmMMBond[bondIter][1] == bondAtom1 ) ) {
997                     
998                     foundBond = true;
999                     break;
1000                 }
1001             }
1002             
1003             if (! foundBond) {
1004                 iout << iERROR << "Error parsing link atom element selection: " 
1005                 << current->data
1006                 << "\n" << endi;
1007                 iout << iERROR << "Atom pair was not found to be a QM-MM bond.\n"
1008                 << endi;
1009                 NAMD_die("Error processing QM information.");
1010             }
1011         }
1012         
1013         strcpy(qmDummyElement[bondIter],linkElement.c_str());
1014         qmDummyElement[bondIter][2] = '\0';
1015         
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]
1020         << "\n" << endi;
1021     }
1022     
1023     
1024     
1025     ///////////////////////////////
1026     /// Bond Schemes. Prepares for treatment of QM-MM bonds in ComputeQM.C
1027     
1028     
1029     int32 **bondsWithAtomLocal = NULL ;
1030     int32 *byAtomSizeLocal = NULL;
1031     ObjectArena <int32 >* tmpArenaLocal = NULL;
1032     if (qmNumBonds > 0) {
1033         
1034         bondsWithAtomLocal = new int32 *[numAtoms];
1035         byAtomSizeLocal = new int32[numAtoms];
1036         tmpArenaLocal = new ObjectArena<int32>;
1037         
1038         //  Build the bond lists
1039        for (int i=0; i<numAtoms; i++)
1040        {
1041          byAtomSizeLocal[i] = 0;
1042        }
1043        for (int i=0; i<numRealBonds; i++)
1044        {
1045          byAtomSizeLocal[bonds[i].atom1]++;
1046          byAtomSizeLocal[bonds[i].atom2]++;
1047        }
1048        for (int i=0; i<numAtoms; i++)
1049        {
1050          bondsWithAtomLocal[i] = tmpArenaLocal->getNewArray(byAtomSizeLocal[i]+1);
1051          bondsWithAtomLocal[i][byAtomSizeLocal[i]] = -1;
1052          byAtomSizeLocal[i] = 0;
1053        }
1054        for (int i=0; i<numRealBonds; i++)
1055        {
1056          int a1 = bonds[i].atom1;
1057          int a2 = bonds[i].atom2;
1058          bondsWithAtomLocal[a1][byAtomSizeLocal[a1]++] = i;
1059          bondsWithAtomLocal[a2][byAtomSizeLocal[a2]++] = i;
1060        }
1061     }
1062     
1063     // In this loops we try to find other bonds in which the MM atoms from
1064     // QM-MM bonds may be involved. The other MM atoms (which we will call M2 and M3)
1065     // will be involved in charge manipulation. See ComputeQM.C for details.
1066     for (int qmBndIt = 0; qmBndIt < qmNumBonds; qmBndIt++) {
1067         
1068         // The charge targets are accumulated in a temporary vector and then 
1069         // transfered to an array that will be transmited to the ComputeQMMgr object.
1070         std::vector<int> chrgTrgt ;
1071         
1072         int MM1 = qmMMBond[qmBndIt][0], MM2, MM2BondIndx, MM3, MM3BondIndx;
1073         
1074         switch (simParams->qmBondScheme) {
1075             
1076             case QMSCHEMERCD:
1077             
1078             case QMSCHEMECS:
1079             {
1080                 // Selects ALL MM2 atoms.
1081                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
1082                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1083                     
1084                     // Checks which of the atoms in the bond structure is the
1085                     // MM2 atom.
1086                     if (bonds[MM2BondIndx].atom1 == MM1)
1087                         MM2 = bonds[MM2BondIndx].atom2;
1088                     else
1089                         MM2 = bonds[MM2BondIndx].atom1;
1090                     
1091                     // In case the bonded atom is a QM atom,
1092                     // skips the index.
1093                     if (qmAtomGroup[MM2] > 0)
1094                         continue;
1095                     
1096                     chrgTrgt.push_back(MM2);
1097                 }
1098                 
1099             } break;
1100             
1101             case QMSCHEMEZ3:
1102             {
1103                 // Selects all MM3 atoms
1104                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
1105                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1106                     
1107                     // Checks which of the atoms in the bond structure is the
1108                     // MM2 atom.
1109                     if (bonds[MM2BondIndx].atom1 == MM1)
1110                         MM2 = bonds[MM2BondIndx].atom2;
1111                     else
1112                         MM2 = bonds[MM2BondIndx].atom1;
1113                     
1114                     // In case the bonded atom is a QM atom,
1115                     // skips the index.
1116                     if (qmAtomGroup[MM2] > 0)
1117                         continue;
1118                     
1119                     for (int i=0; i<byAtomSizeLocal[MM2]; i++) {
1120                         MM3BondIndx = bondsWithAtomLocal[MM2][i] ;
1121                         
1122                         // Checks which of the atoms in the bond structure is the
1123                         // MM3 atom.
1124                         if (bonds[MM3BondIndx].atom1 == MM2)
1125                             MM3 = bonds[MM3BondIndx].atom2;
1126                         else
1127                             MM3 = bonds[MM3BondIndx].atom1;
1128                         
1129                         // In case the bonded atom is a QM atom,
1130                         // skips the index.
1131                         // We also keep the search from going back to MM1.
1132                         if (qmAtomGroup[MM3] > 0 || MM3 == MM1)
1133                             continue;
1134                         
1135                         chrgTrgt.push_back(MM3);
1136                     }
1137                     
1138                 }
1139                 
1140             };
1141             
1142             case QMSCHEMEZ2:
1143             {
1144                 // Selects all MM2 atoms
1145                 for (int i=0; i<byAtomSizeLocal[MM1]; i++) {
1146                     MM2BondIndx = bondsWithAtomLocal[MM1][i] ;
1147                     
1148                     // Checks which of the atoms in the bond structure is the
1149                     // MM2 atom.
1150                     if (bonds[MM2BondIndx].atom1 == MM1)
1151                         MM2 = bonds[MM2BondIndx].atom2;
1152                     else
1153                         MM2 = bonds[MM2BondIndx].atom1;
1154                     
1155                     // In case the bonded atom is a QM atom,
1156                     // skips the index.
1157                     if (qmAtomGroup[MM2] > 0)
1158                         continue;
1159                     
1160                     chrgTrgt.push_back(MM2);
1161                 }
1162                 
1163             };
1164             
1165             case QMSCHEMEZ1:
1166             {
1167                 // Selects all MM1 atoms
1168                 chrgTrgt.push_back(MM1);
1169             } break;
1170         }
1171         
1172         
1173         qmMMChargeTarget[qmBndIt] = new int[chrgTrgt.size()] ;
1174         qmMMNumTargs[qmBndIt] =  chrgTrgt.size();
1175         
1176         DebugM(3, "MM-QM bond " << qmBndIt << "; MM atom " 
1177         << qmMMBond[qmBndIt][0] << " conections: \n" );
1178         
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" );
1182         }
1183         
1184         chrgTrgt.clear();
1185     }
1186     
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;
1193     
1194     
1195     ///////////////////////////////
1196     /// Live Solvent Selection
1197     
1198     
1199     if(simParams->qmLSSOn) {
1200         
1201         std::map<Real,int> grpLSSSize ;
1202         std::map<Real,int>::iterator itGrpSize;
1203         
1204         qmLSSTotalNumAtms = 0;
1205         qmLSSResidueSize = 0;
1206         
1207         if (simParams->qmLSSFreq == 0)
1208             qmLSSFreq = simParams->stepsPerCycle ;
1209         else
1210             qmLSSFreq = simParams->qmLSSFreq;
1211             
1212         #ifdef DEBUG_QM
1213         int resSize = -1;
1214         #endif
1215         
1216         std::map<Real, int> grpNumLSSRes;
1217         std::map<Real, int>::iterator itGrpNumRes;
1218         
1219         for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
1220             
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" 
1224                 << endi;
1225                 NAMD_die("Error parsing data for LSS.");
1226             }
1227             
1228             qmLSSTotalNumAtms += it->size;
1229             
1230             #ifdef DEBUG_QM
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" 
1235                 << endi;
1236                 NAMD_die("Error parsing data for LSS.");
1237             }
1238                 
1239 //             DebugM(3,"Residue " << it->resID << ": " << it->segName
1240 //             << " - from " << it->begAtmID << " with size "
1241 //             << it->size << " (QM ID: " << it->qmGrpID
1242 //             << ") has " << it->atmIDs.size() << " atoms: \n" ) ;
1243 //             for (int i=0; i<it->atmIDs.size(); i++) 
1244 //                 DebugM(3, it->atmIDs[i] << "\n" );
1245             #endif
1246             
1247             // Calculating total number of atoms per group
1248             itGrpSize = grpLSSSize.find(it->qmGrpID) ;
1249             if (itGrpSize != grpLSSSize.end())
1250                 itGrpSize->second += it->size;
1251             else
1252                 grpLSSSize.insert(std::pair<Real,int>(it->qmGrpID, it->size));
1253             
1254             // Calculating total number of solvent residues per group
1255             itGrpNumRes = grpNumLSSRes.find(it->qmGrpID) ;
1256             if (itGrpNumRes != grpNumLSSRes.end())
1257                 itGrpNumRes->second += 1;
1258             else
1259                 grpNumLSSRes.insert(std::pair<Real,int>(it->qmGrpID, 1));
1260         }
1261         
1262         qmLSSResidueSize = lssGrpRes.begin()->size ;
1263         
1264         qmLSSSize = new int[qmNumGrps];
1265         
1266         qmLSSIdxs = new int[qmLSSTotalNumAtms];
1267         int lssAtmIndx = 0;
1268         
1269         switch (simParams->qmLSSMode) {
1270         
1271         case QMLSSMODECOM:
1272         {
1273             
1274             qmLSSRefSize = new int[qmNumGrps];
1275             
1276             qmLSSMass = new Mass[qmLSSTotalNumAtms];
1277             
1278             qmLSSRefIDs = new int[totalNumRefAtms];
1279             qmLSSRefMass = new Mass[totalNumRefAtms];
1280             int lssRefIndx = 0;
1281             
1282             for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
1283                 
1284                 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
1285                 
1286                 if (itGrpSize != grpNumLSSRes.end())
1287                     qmLSSSize[grpIndx] =  itGrpSize->second;
1288                 else
1289                     qmLSSSize[grpIndx] =  0;
1290                 
1291                 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
1292                     
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;
1297                             lssAtmIndx++;
1298                         }
1299                     }
1300                 }
1301                 
1302                 DebugM(4, "QM group " << qmGrpID[grpIndx]
1303                 << " has " << qmLSSSize[grpIndx] << " solvent molecules, "
1304                 << "totalling " << grpLSSSize[qmGrpID[grpIndx]]
1305                 << " atoms (check " << lssAtmIndx << ").\n" );
1306                 
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;
1311                     lssRefIndx++;
1312                 }
1313                 
1314                 DebugM(4, "QM group " << qmGrpID[grpIndx]
1315                 << " has " << qmLSSRefSize[grpIndx] << " non-solvent atoms for LSS.\n" );
1316             }
1317             
1318         } break ;
1319         
1320         case QMLSSMODEDIST:
1321         {
1322             for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
1323                 
1324                 itGrpSize = grpNumLSSRes.find(qmGrpID[grpIndx]) ;
1325                 
1326                 if (itGrpSize != grpNumLSSRes.end())
1327                     qmLSSSize[grpIndx] =  itGrpSize->second;
1328                 else
1329                     qmLSSSize[grpIndx] =  0;
1330                 
1331                 for( auto it = lssGrpRes.begin(); it != lssGrpRes.end();it++ ) {
1332                     
1333                     if (it->qmGrpID == qmGrpID[grpIndx]) {
1334                         for (int i=0; i<it->atmIDs.size(); i++) {
1335                             qmLSSIdxs[lssAtmIndx] = it->atmIDs[i];
1336                             lssAtmIndx++;
1337                         }
1338                     }
1339                 }
1340                 
1341                 DebugM(4, "QM group " << qmGrpID[grpIndx]
1342                 << " has " << qmLSSSize[grpIndx] << " solvent molecules, "
1343                 << "totalling " << grpLSSSize[qmGrpID[grpIndx]]
1344                 << " atoms (check " << lssAtmIndx << ").\n" );
1345             }
1346             
1347         } break ;
1348             
1349         }
1350     }
1351     
1352     
1353     ///////////////////////////////
1354     /// Custom Point Charge selection
1355     
1356     
1357     PDB *customPCPDB;
1358     
1359     // In case we have a custom and fixed set of point charges for each QM group,
1360     // we process the files containing information.
1361     current = NULL;
1362     if (simParams->qmCustomPCSel) {
1363         current = cfgList->find("QMCustomPCFile");
1364     }
1365     
1366     std::map<Real,std::vector<int> > qmPCVecMap ;
1367     
1368     int numParsedPBDs = 0;
1369     for ( ; current != NULL; current = current->next ) {
1370         
1371         iout << iINFO << "Parsing QM Custom PC file " << current->data << "\n" << endi;
1372         customPCPDB = new PDB(current->data);
1373         
1374         if (customPCPDB->num_atoms() != numAtoms)
1375            NAMD_die("Number of atoms in QM Custom PC PDB file doesn't match coordinate PDB");
1376         
1377         std::vector< int > currPCSel ;
1378         Real currQMID = 0 ;
1379         int currGrpSize = 0 ;
1380         
1381         for (int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
1382             
1383             BigReal beta = customPCPDB->atom(atmInd)->temperaturefactor() ;
1384             BigReal occ = customPCPDB->atom(atmInd)->occupancy() ;
1385             
1386             if ( beta != 0 && occ != 0)
1387                 NAMD_die("An atom cannot be marked as QM and poitn charge simultaneously!");
1388             
1389             // If this is not a QM atom and 
1390             if (occ != 0) {
1391                 currPCSel.push_back(atmInd) ;
1392             }
1393             
1394             if (beta != 0) {
1395                 if (pdbP->atom(atmInd)->temperaturefactor() != beta)
1396                     NAMD_die("QM Group selection is different in reference PDB and Custom-PC PDB!");
1397                 
1398                 if (currQMID == 0) {
1399                     // If this is the first QM atom we find, keep the QM Group ID.
1400                     currQMID = beta;
1401                 }
1402                 else {
1403                     // For all other atoms, check if it is the same group. It must be!!
1404                     if (currQMID != beta)
1405                         NAMD_die("Found two different QM group IDs in this file!");
1406                 }
1407                 
1408                 currGrpSize++;
1409             }
1410             
1411         }
1412         
1413         if (currGrpSize != qmGrpSizeVec[ qmGrpIDMap[currQMID] ])
1414             NAMD_die("Reference PDB and Custom-PC PDB have different QM group sizes!") ;
1415         
1416         qmPCVecMap.insert(std::pair<Real,std::vector<int> >(
1417             currQMID, currPCSel ));
1418         
1419         numParsedPBDs++;
1420         delete customPCPDB;
1421     }
1422     
1423     delete pdbP;
1424     
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;
1428     }
1429     
1430     // Initializes an array with the number of Custom Point Charges per 
1431     // QM group.
1432     qmCustPCSizes = new int[qmNumGrps];
1433     for (int i=0; i<qmNumGrps; i++)
1434         qmCustPCSizes[i] = 0;
1435     
1436     qmTotCustPCs = 0;
1437     
1438     // Stores the size of each Custom PC vector in the array.
1439     // We may not have PCs for all QM groups.
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();
1444     }
1445     
1446     qmCustomPCIdxs = new int[qmTotCustPCs];
1447     
1448     if (simParams->qmCustomPCSel) {
1449         
1450         int auxIter = 0;
1451         for (int grpIndx=0; grpIndx<qmNumGrps; grpIndx++) {
1452             
1453             Real currQMID = qmGrpID[grpIndx];
1454             
1455             iout << iINFO << "Loading " << qmPCVecMap[currQMID].size()
1456             << " custom point charges to QM Group " << grpIndx
1457             << " (ID: " << currQMID << ")\n" << endi;
1458             
1459             for (int pcIndx=0; pcIndx<qmPCVecMap[currQMID].size(); pcIndx++) {
1460                 qmCustomPCIdxs[auxIter] = qmPCVecMap[currQMID][pcIndx] ;
1461                 auxIter++;
1462             }
1463         }
1464     } 
1465     
1466     // Conditional SMD option
1467     qmcSMD = simParams->qmCSMD;
1468     if (qmcSMD) {
1469         read_qm_csdm_file(qmGrpIDMap);
1470     }
1471     
1472     ///////////////////////////////
1473     /// Topology preparation
1474     
1475     if (simParams->qmElecEmbed) {
1476         // Modifies Atom charges for Electrostatic Embedding.
1477         // QM atoms cannot have charges in the standard location, to keep
1478         // NAMD from calculating electrostatic interactions between QM and MM atoms.
1479         // We handle electrostatics ourselves in ComputeQM.C and in special
1480         // modifications for PME.
1481         for (int i=0; i<qmNumQMAtoms; i++) {
1482             qmAtmChrg[i] = atoms[qmAtmIndx[i]].charge;
1483             atoms[qmAtmIndx[i]].charge = 0;
1484         }
1485     }
1486     
1487     
1488     if ( simParams->extraBondsOn) {
1489         // Lifted from Molecule::build_extra_bonds
1490         
1491         StringList *file = cfgList->find("extraBondsFile");
1492         
1493         char err_msg[512];
1494         int a1,a2; float k, ref;
1495         
1496         for ( ; file; file = file->next ) {  // loop over files
1497             FILE *f = fopen(file->data,"r");
1498 //             if ( ! f ) {
1499 //               sprintf(err_msg, "UNABLE TO OPEN EXTRA BONDS FILE %s", file->data);
1500 //               NAMD_err(err_msg);
1501 //             } else {
1502 //               iout << iINFO << "READING EXTRA BONDS FILE " << file->data <<"\n"<<endi;
1503 //             }
1504             
1505             while ( 1 ) {
1506               char buffer[512];
1507               int ret_code;
1508               do {
1509                 ret_code = NAMD_read_line(f, buffer);
1510               } while ( (ret_code==0) && (NAMD_blank_string(buffer)) );
1511               if (ret_code!=0) break;
1512
1513               char type[512];
1514               sscanf(buffer,"%s",type);
1515               
1516               int badline = 0;
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;
1520                 
1521                 // If an extra bond is defined between QM atoms, we make
1522                 // note so that it wont be deleted when we delete bonded 
1523                 // interactions between QM atoms.
1524                 if( qmAtomGroup[a1] > 0 && qmAtomGroup[a2]) {
1525                     Bond tmp;
1526                     tmp.bond_type = 0;
1527                     tmp.atom1 = a1;  tmp.atom2 = a2;
1528                     qmExtraBonds.add(tmp);
1529                 }
1530                 
1531                 
1532               }else if ( ! strncasecmp(type,"#",1) ) {
1533                 continue;  // comment
1534               }
1535               
1536             }
1537             fclose(f);
1538         }
1539     }
1540     
1541     return;
1542     
1543 #endif // MEM_OPT_VERSION
1544 }
1545
1546 // Adapted from Molecule::delete_alch_bonded
1547 void Molecule::delete_qm_bonded(void)  {
1548
1549 #ifdef MEM_OPT_VERSION
1550     NAMD_die("QMMM interface is not supported in memory optimized builds.");
1551 #else
1552     
1553     DebugM(3,"Cleaning QM bonds, angles, dihedrals, impropers and crossterms.\n")
1554     
1555     // Bonds
1556     Bond *nonQMBonds;
1557     nonQMBonds = new Bond[numBonds] ;
1558     int nonQMBondCount = 0;
1559     qmDroppedBonds = 0;
1560     for (int i = 0; i < numBonds; i++) {
1561         
1562         int part1 = qmAtomGroup[bonds[i].atom1];
1563         int part2 = qmAtomGroup[bonds[i].atom2];
1564         if (part1 > 0 && part2 > 0 ) {
1565             
1566             // If the user defined extra bonds, we check if the QM-QM bond is an "extra"
1567             // bond, and do not delete it.
1568             if (simParams->extraBondsOn) {
1569 //                 std::cout << "Checking Bond: " << bonds[i].atom1 << "," << bonds[i].atom2 << std::endl ;
1570                 
1571                 for (int ebi=0; ebi < qmExtraBonds.size() ; ebi++) {
1572                     
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) ) {
1577 //                         std::cout << "This is an extra bond! We will keep it." << std::endl;
1578                         nonQMBonds[nonQMBondCount++] = bonds[i];
1579                         break;
1580                     }
1581                 }
1582             }
1583             
1584             qmDroppedBonds++;
1585         } else {
1586             // Just a simple sanity check.
1587             // If there are no QM bonds defined by the user but we find a
1588             // bond between a QM and an MM atom, error out.
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.");
1593             }
1594             nonQMBonds[nonQMBondCount++] = bonds[i];
1595         }
1596     }
1597     numBonds = nonQMBondCount;
1598     delete [] bonds;
1599     bonds = new Bond[numBonds];
1600     for (int i = 0; i < nonQMBondCount; i++) {
1601         bonds[i]=nonQMBonds[i];
1602     }
1603     delete [] nonQMBonds;
1604   numRealBonds = numBonds;
1605     
1606   // Angles
1607   Angle *nonQMAngles;
1608   nonQMAngles = new Angle[numAngles];
1609   int nonQMAngleCount = 0;
1610   qmDroppedAngles = 0;
1611   for (int i = 0; i < numAngles; i++) {
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) {
1616       //CkPrintf("-----ANGLE ATOMS %i %i %i partitions %i %i %i\n",angles[i].atom1, angles[i].atom2, angles[i].atom3, part1, part2, part3);
1617       qmDroppedAngles++;
1618     }
1619     else {
1620       nonQMAngles[nonQMAngleCount++] = angles[i];
1621     }
1622   }
1623   numAngles = nonQMAngleCount;
1624   delete [] angles;
1625   angles = new Angle[numAngles];
1626   for (int i = 0; i < nonQMAngleCount; i++) {
1627     angles[i]=nonQMAngles[i];
1628   }
1629   delete [] nonQMAngles;
1630
1631
1632   // Dihedrals
1633   Dihedral *nonQMDihedrals;
1634   nonQMDihedrals = new Dihedral[numDihedrals];
1635   int nonQMDihedralCount = 0;
1636   qmDroppedDihedrals = 0;
1637   for (int i = 0; i < numDihedrals; i++) {
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) {
1643       //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);
1644       qmDroppedDihedrals++;
1645     }
1646     else {
1647       nonQMDihedrals[nonQMDihedralCount++] = dihedrals[i];
1648     }
1649   }
1650   numDihedrals = nonQMDihedralCount;
1651   delete [] dihedrals;
1652   dihedrals = new Dihedral[numDihedrals];
1653   for (int i = 0; i < numDihedrals; i++) {
1654     dihedrals[i]=nonQMDihedrals[i];
1655   }
1656   delete [] nonQMDihedrals;
1657
1658   // Impropers
1659   Improper *nonQMImpropers;
1660   nonQMImpropers = new Improper[numImpropers];
1661   int nonQMImproperCount = 0;
1662   qmDroppedImpropers = 0;
1663   for (int i = 0; i < numImpropers; i++) {
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) {
1669       //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);
1670       qmDroppedImpropers++;
1671     }
1672     else {
1673       nonQMImpropers[nonQMImproperCount++] = impropers[i];
1674     }
1675   }
1676   numImpropers = nonQMImproperCount;
1677   delete [] impropers;
1678   impropers = new Improper[numImpropers];
1679   for (int i = 0; i < numImpropers; i++) {
1680     impropers[i]=nonQMImpropers[i];
1681   }
1682   delete [] nonQMImpropers;
1683   
1684   // Crossterms 
1685   Crossterm *nonQMCrossterms;
1686   nonQMCrossterms = new Crossterm[numCrossterms];
1687   int nonQMCrosstermCount = 0;
1688   qmDroppedCrossterms = 0 ;
1689   for (int i = 0; i < numCrossterms; i++) {
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) {
1700       //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);
1701       qmDroppedCrossterms++;
1702     }
1703     else {
1704       nonQMCrossterms[nonQMCrosstermCount++] = crossterms[i];
1705     }
1706   }
1707   numCrossterms = nonQMCrosstermCount;
1708   delete [] crossterms;
1709   crossterms = new Crossterm[numCrossterms] ;
1710   for (int i=0; i < numCrossterms; i++){
1711       crossterms[i] = nonQMCrossterms[i] ;
1712   }
1713   delete [] nonQMCrossterms ;
1714   
1715   if (!CkMyPe()) {
1716       iout << iINFO << "The QM region will remove " << qmDroppedBonds << " bonds, " << 
1717           qmDroppedAngles << " angles, " << qmDroppedDihedrals << " dihedrals, "
1718           << qmDroppedImpropers << " impropers and " << qmDroppedCrossterms 
1719           << " crossterms.\n" << endi ;
1720   }
1721   
1722 #endif
1723 }
1724
1725 typedef std::pair<int,int> cSMDPair ;
1726
1727 void Molecule::read_qm_csdm_file(std::map<Real,int> &qmGrpIDMap)  {
1728     
1729     std::ifstream cSMDInputFile ;
1730     std::string Line("") ;
1731     
1732     iout << iINFO << "Reading QM cSMD configuration file: " << 
1733         simParams->qmCSMDFile << "\n" << endi ;
1734     
1735     cSMDInputFile.open( simParams->qmCSMDFile ) ;
1736     
1737     if (! cSMDInputFile.good())
1738     {
1739         NAMD_die("Configuration file for QM cSMD could not be read!") ;
1740     }
1741     
1742     
1743     // Index of Conditional SMD guides per group
1744     ResizeArray<ResizeArray<int> > cSMDindexV;
1745     cSMDindexV.resize(qmNumGrps);
1746     // Atom indices for Origin and Target of cSMD
1747     ResizeArray<cSMDPair> cSMDpairsV;
1748     // Spring constants for cSMD
1749     ResizeArray<Real> cSMDKsV;
1750     // Speed of movement of guide particles for cSMD.
1751     ResizeArray<Real> cSMDVelsV;
1752     // Distance cutoff for guide particles for cSMD.
1753     ResizeArray<Real> cSMDcoffsV;
1754     
1755     while( ! cSMDInputFile.eof() )
1756     {
1757         getline(cSMDInputFile, Line) ;
1758         
1759         if ( ! Line.length() )
1760             continue;
1761         
1762         if (Line.substr(0,1) == std::string("#") )
1763             continue;
1764         
1765         auto strVec = split( Line, " ");
1766         
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.");
1771         }
1772         
1773         std::stringstream storConv ;
1774         Real convData ;
1775         
1776         cSMDpairsV.add( cSMDPair(0,0) );
1777         cSMDKsV.add(0);
1778         cSMDVelsV.add(0);
1779         cSMDcoffsV.add(0);
1780         
1781         for (int i=0; i < 5; i++ ) {
1782             storConv.clear() ;
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.");
1789             }
1790             
1791             switch (i) {
1792             case 0:
1793                 cSMDpairsV[cSMDnumInst].first = convData ;
1794                 break;
1795             case 1:
1796                 cSMDpairsV[cSMDnumInst].second = convData ;
1797                 break;
1798             case 2:
1799                 cSMDKsV[cSMDnumInst] = convData ;
1800                 break;
1801             case 3:
1802                 cSMDVelsV[cSMDnumInst] = convData ;
1803                 break;
1804             case 4:
1805                 cSMDcoffsV[cSMDnumInst] = convData ;
1806                 break;
1807             }
1808         }
1809         
1810         // Sanity check
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.") ;
1815         }
1816         
1817         // Sanity check
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.") ;
1823         }
1824         
1825         // Sanity check
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.") ;
1830         }
1831         
1832         int grpIndx = qmGrpIDMap[qmAtomGroup[cSMDpairsV[cSMDnumInst].first]] ;
1833         
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 ;
1838         
1839         cSMDindexV[grpIndx].add(cSMDnumInst);
1840         
1841         cSMDnumInst++;
1842     }
1843     
1844     cSMDindex = new int*[qmNumGrps];
1845     cSMDindxLen = new int[qmNumGrps];
1846     
1847     for (int i=0; i<qmNumGrps; i++) {
1848             cSMDindex[i] = new int[cSMDindexV[i].size()];
1849             cSMDindxLen[i] = cSMDindexV[i].size();
1850     
1851             for (int j=0; j<cSMDindxLen[i]; j++)
1852                     cSMDindex[i][j] = cSMDindexV[i][j] ;
1853     }
1854     
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];
1861
1862     for (int i=0; i<cSMDnumInst; i++) {
1863             cSMDpairs[i][0] = cSMDpairsV[i].first;
1864             cSMDpairs[i][1] = cSMDpairsV[i].second;
1865     
1866             cSMDKs[i] = cSMDKsV[i];
1867             cSMDVels[i] = cSMDVelsV[i];
1868             cSMDcoffs[i] = cSMDcoffsV[i];
1869     }
1870     
1871 }
1872
1873
1874
1875
1876