NAMD
MoleculeQM.C
Go to the documentation of this file.
1 
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];
50  std::vector<int> atmIDs ;
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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
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;
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 
char qmLSSResname[5]
int get_residue_size(const char *segid, int resid) const
Definition: Molecule.C:144
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Type * getNewArray(int n)
Definition: ObjectArena.h:49
std::map< Real, refSelStrVec > refSelStrMap
Definition: MoleculeQM.C:106
#define QMLSSMODECOM
Definition: Molecule.h:126
int begAtmID
Definition: MoleculeQM.C:49
Definition: PDB.h:35
int numBonds
Definition: Molecule.h:559
#define QMSCHEMEZ2
Definition: Molecule.h:139
int NAMD_read_line(FILE *fd, char *buf, int bufsize)
Definition: strlib.C:38
void assign_vdw_index(const char *, Atom *)
Definition: Parameters.C:3470
std::vector< int > atmIDs
Definition: MoleculeQM.C:50
short int32
Definition: dumpdcd.c:24
#define QMSCHEMERCD
Definition: Molecule.h:137
BigReal temperaturefactor(void)
Definition: PDBData.C:450
float Real
Definition: common.h:109
#define DebugM(x, y)
Definition: Debug.h:59
char qmCSMDFile[NAMD_FILENAME_BUFFER_SIZE]
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
int numRealBonds
Definition: Molecule.h:558
int32 atom1
Definition: structures.h:81
refSelStr(std::string newSegID, int newResid)
Definition: MoleculeQM.C:99
#define QMFormatMOPAC
Definition: Molecule.h:129
#define iout
Definition: InfoStream.h:87
int num_atoms(void)
Definition: PDB.C:323
#define QMFormatORCA
Definition: Molecule.h:128
int insert(const Elem &elem)
#define QMSCHEMEZ3
Definition: Molecule.h:140
int32 atom1
Definition: structures.h:48
#define QMSCHEMEZ1
Definition: Molecule.h:138
#define QMSCHEMECS
Definition: Molecule.h:136
const char * residuename(void)
Definition: PDBData.C:399
PDBAtom * atom(int place)
Definition: PDB.C:394
int NAMD_blank_string(char *str)
Definition: strlib.C:222
iterator end(void)
Definition: ResizeArray.h:37
Real qmGrpID
Definition: MoleculeQM.C:51
int32 atom1
Definition: structures.h:55
int Bool
Definition: common.h:133
int residueseq(void)
Definition: PDBData.C:411
std::string segid
Definition: MoleculeQM.C:95
std::pair< Real, refSelStrVec > refSelStrPair
Definition: MoleculeQM.C:107
int lookup(const char *segid, int resid, int *begin, int *end) const
Definition: Molecule.C:76
void delete_qm_bonded()
Definition: MoleculeQM.C:1547
void prepare_qm(const char *pdbFileName, Parameters *params, ConfigList *cfgList)
Definition: MoleculeQM.C:109
int numAngles
Definition: Molecule.h:560
int numAtoms
Definition: Molecule.h:556
int numCrossterms
Definition: Molecule.h:567
void NAMD_die(const char *err_msg)
Definition: common.C:83
bool operator<(const qmSolvData &ref)
Definition: MoleculeQM.C:69
int resid
Definition: MoleculeQM.C:96
std::vector< std::string > split(const std::string &text, std::string delimiter)
Definition: MoleculeQM.C:73
int numDihedrals
Definition: Molecule.h:561
int numImpropers
Definition: Molecule.h:566
int add(const Elem &elem)
Definition: ResizeArray.h:97
int32 atom1
Definition: structures.h:72
Bool qmMOPACAddConfigChrg
std::vector< refSelStr > refSelStrVec
Definition: MoleculeQM.C:105
StringList * next
Definition: ConfigList.h:49
void resize(int i)
Definition: ResizeArray.h:84
char * data
Definition: ConfigList.h:48
char segName[5]
Definition: MoleculeQM.C:48
bool operator==(const qmSolvData &ref)
Definition: MoleculeQM.C:70
StringList * find(const char *name) const
Definition: ConfigList.C:341
Elem * find(const Elem &elem)
std::pair< Origin, Target > cSMDPair
Definition: ComputeQM.h:112
float Mass
Definition: ComputeGBIS.inl:20
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
BigReal occupancy(void)
Definition: PDBData.C:444
int size(void) const
Definition: ResizeArray.h:127
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
qmSolvData(int newResID, int newBegAtm, int newSize, const char *newSegName, Real newQMID)
Definition: MoleculeQM.C:59
#define QMLSSMODEDIST
Definition: Molecule.h:125
Index bond_type
Definition: structures.h:50
const char * element(void)
Definition: PDBData.C:470
int32 atom1
Definition: structures.h:63
int32 atom2
Definition: structures.h:49
double BigReal
Definition: common.h:114
char qmColumn[16]
const char * segmentname(void)
Definition: PDBData.C:464
qmSolvData(int newResID, int newBegAtm, int newSize)
Definition: MoleculeQM.C:54
iterator begin(void)
Definition: ResizeArray.h:36