NAMD
ComputeQM.C
Go to the documentation of this file.
1 
7 #include "dcdlib.h"
8 
9 // #define DEBUGM
10 // #define DEBUG_QM
11 
12 #ifdef DEBUG_QM
13  #define DEBUGM
14 #endif
15 
16 #define MIN_DEBUG_LEVEL 1
17 #include "Debug.h"
18 
19 #include "InfoStream.h"
20 #include "Node.h"
21 #include "PatchMap.h"
22 #include "PatchMap.inl"
23 #include "AtomMap.h"
24 #include "ComputeQM.h"
25 #include "ComputeQMMgr.decl.h"
26 #include "PatchMgr.h"
27 #include "Molecule.h"
28 #include "ReductionMgr.h"
29 #include "ComputeMgr.h"
30 #include "ComputeMgr.decl.h"
31 #include "SimParameters.h"
32 #include "WorkDistrib.h"
33 #include "varsizemsg.h"
34 #include <stdlib.h>
35 #include <stdio.h>
36 #include <errno.h>
37 #include <algorithm>
38 
39 #include "ComputePme.h"
40 #include "ComputePmeMgr.decl.h"
41 
42 #include <fstream>
43 #include <iomanip>
44 
45 #if defined(WIN32) && !defined(__CYGWIN__)
46 #include <direct.h>
47 #define mkdir(X,Y) _mkdir(X)
48 #define S_ISDIR(X) ((X) & S_IFDIR)
49 #endif
50 
51 #ifndef SQRT_PI
52 #define SQRT_PI 1.7724538509055160273 /* mathematica 15 digits*/
53 #endif
54 
55 #define QMLENTYPE 1
56 #define QMRATIOTYPE 2
57 
58 #define QMPCSCHEMENONE 1
59 #define QMPCSCHEMEROUND 2
60 #define QMPCSCHEMEZERO 3
61 
62 #define QMPCSCALESHIFT 1
63 #define QMPCSCALESWITCH 2
64 
65 struct ComputeQMAtom {
67  float charge;
68  int id;
70  int homeIndx;
71  int vdwType;
72  ComputeQMAtom() : position(0), charge(0), id(-1), qmGrpID(-1),
73  homeIndx(-1), vdwType(0) {}
74  ComputeQMAtom(Position posInit, float chrgInit, int idInit,
75  Real qmInit, int hiInit, int newvdwType) {
76  position = posInit;
77  charge = chrgInit;
78  id = idInit;
79  qmGrpID = qmInit;
80  homeIndx = hiInit;
81  vdwType = newvdwType;
82  }
84  position = ref.position;
85  charge = ref.charge;
86  id = ref.id;
87  qmGrpID = ref.qmGrpID;
88  homeIndx = ref.homeIndx;
89  vdwType = ref.vdwType;
90  }
91 };
92 
93 // sort using a custom function
95 {
96  return a.id < b.id;
97 }
98 
99 #define PCMODEUPDATESEL 1
100 #define PCMODEUPDATEPOS 2
101 #define PCMODECUSTOMSEL 3
102 
103 class QMCoordMsg : public CMessage_QMCoordMsg {
104 public:
106  int numAtoms;
107  int timestep;
111  int *pcIndxs;
112 };
113 
114 struct pntChrgDist {
115  int index;
117  pntChrgDist() : index(-1), dist(0) {};
118  pntChrgDist(int newIndx, Real newDist) {
119  index = newIndx;
120  dist = newDist;
121  }
122  int operator <(const pntChrgDist& v) {return dist < v.dist;}
123 };
124 
127  float charge;
128  int id;
130  int homeIndx;
133  int vdwType;
134  ComputeQMPntChrg() : position(0), charge(0), id(-1), qmGrpID(-1),
135  homeIndx(-1), dist(0), mass(0), vdwType(0) {}
136  ComputeQMPntChrg(Position posInit, float chrgInit, int idInit,
137  Real qmInit, int hiInit, Real newDist,
138  Mass newM, int newvdwType) {
139  position = posInit;
140  charge = chrgInit;
141  id = idInit;
142  qmGrpID = qmInit;
143  homeIndx = hiInit;
144  dist = newDist;
145  mass = newM;
146  vdwType = newvdwType;
147  }
149  position = ref.position;
150  charge = ref.charge;
151  id = ref.id;
152  qmGrpID = ref.qmGrpID;
153  homeIndx = ref.homeIndx;
154  dist = ref.dist;
155  mass = ref.mass;
156  vdwType = ref.vdwType;
157  }
158 
159  bool operator <(const ComputeQMPntChrg& ref) {
160  return (id < ref.id);
161  }
162  bool operator ==(const ComputeQMPntChrg& ref) {
163  return (id == ref.id) ;
164  }
165 };
166 
167 class QMPntChrgMsg : public CMessage_QMPntChrgMsg {
168 public:
170  int numAtoms;
172 };
173 
174 class QMGrpResMsg : public CMessage_QMGrpResMsg {
175 public:
176  int grpIndx;
177  BigReal energyOrig; // Original QM Energy
178  BigReal energyCorr; // Corrected energy due to PME
182 };
183 
184 class QMForceMsg : public CMessage_QMForceMsg {
185 public:
186  bool PMEOn;
193 };
194 
195 # define QMATOMTYPE_DUMMY 0
196 # define QMATOMTYPE_QM 1
197 # define QMPCTYPE_IGNORE 0
198 # define QMPCTYPE_CLASSICAL 1
199 # define QMPCTYPE_EXTRA 2
200 
201 #define QMLSSQMRES 1
202 #define QMLSSCLASSICALRES 2
203 
204 struct idIndxStr {
205  int ID; // Global atom ID
206  int indx; // Atom index in the vector received from other ranks.
207 
209  ID = -1;
210  indx = -1;
211  };
212  idIndxStr(int newID, int newIndx) {
213  ID = newID;
214  indx = newIndx;
215  }
216  idIndxStr(const idIndxStr& ref) {
217  ID = ref.ID ;
218  indx = ref.indx;
219  }
220 
222  ID = ref.ID ;
223  indx = ref.indx;
224  return *this;
225  }
226 
227  bool operator==(const idIndxStr& ref) const {
228  return (ID == ref.ID && indx == ref.indx);
229  }
230  bool operator!=(const idIndxStr& ref) const {
231  return !(*this == ref);
232  }
233 
234  // We sort the atoms based on their global ID so that the same
235  // order is kept when replacing classical and QM residues.
236  // We are assuming here that every time a residue is added to a PDB
237  // file, its atoms are placed in the same order.
238  bool operator<(const idIndxStr& ref) {return ID < ref.ID;}
239 };
240 
241 struct lssDistSort {
242  int type;
244 // std::vector<int> idVec;
245 // std::vector<int> indxVec;
247 
249  idIndx.resize(0);
250  };
251  lssDistSort(int newType, Real newDist) {
252  type = newType;
253  dist = newDist;
254  idIndx.resize(0);
255  }
256  lssDistSort(const lssDistSort& ref) {
257  type = ref.type;
258  dist = ref.dist;
260  for (int i=0; i<ref.idIndx.size(); i++) idIndx.insert(ref.idIndx[i]) ;
261  idIndx.sort();
262  }
263 
265  type = ref.type;
266  dist = ref.dist;
268  for (int i=0; i<ref.idIndx.size(); i++) idIndx.insert(ref.idIndx[i]) ;
269  idIndx.sort();
270 
271  return *this;
272  }
273 
274  bool operator==(const lssDistSort& ref) {
275  bool returnVal = true;
276 
277  if (! (type == ref.type && dist == ref.dist))
278  return false;
279 
280  if (idIndx.size() != ref.idIndx.size())
281  return false;
282 
283  for (int i=0; i<ref.idIndx.size(); i++) {
284  if (idIndx[i] != ref.idIndx[i])
285  return false;
286  }
287 
288  return returnVal;
289  }
290 
291  bool operator<(const lssDistSort& ref) {return dist < ref.dist;}
292 
293 } ;
294 
295 struct QMAtomData {
297  float charge;
298  int id;
299  int bountToIndx; // So as to implement an "exclude 1-2" functionality.
300  int type;
301  char element[3];
303  QMAtomData() : position(0), charge(0), id(-1), bountToIndx(-1),
304  type(-1), dist(0) {}
305  QMAtomData(Position posInit, float chrgInit, int idInit,
306  int bountToIndxInit, int newType,
307  char *elementInit, Real newDist) {
308  position = posInit;
309  charge = chrgInit;
310  id = idInit;
311  bountToIndx = bountToIndxInit;
312  type = newType;
313  strncpy(element,elementInit,3);
314  dist = newDist;
315  }
316 };
317 
318 class QMGrpCalcMsg: public CMessage_QMGrpCalcMsg {
319 public:
320  int timestep;
321  int grpIndx;
322  int peIter;
324  int numAllAtoms ; // Including dummy atoms.
326  int numAllPntChrgs; // Inlcuding point charges created to handle QM-MM bonds.
329  bool secProcOn ;
330  bool prepProcOn ;
331  bool PMEOn;
332  bool switching ;
335  int pcScheme ;
338  char baseDir[256], execPath[256], secProc[256], prepProc[256];
340  char *configLines;
341 };
342 
343 struct dummyData {
344  // Position of dummy atom.
346  // Indicates to which QM atom is this dummy atom bound.
347  // This is indexed in the context of the number of QM atoms in a QM group.
348  int qmInd ;
349  // Indicates the index of this bond, used to get the element of the dummy atom.
350  int bondIndx;
351  dummyData(Position newPos, int newQmInd, int newIndx) {
352  pos = newPos;
353  qmInd = newQmInd;
354  bondIndx = newIndx;
355  }
356 } ;
357 
358 struct LSSDataStr {
359  int resIndx ;
361  LSSDataStr(int newResIndx, Mass newMass) {
362  resIndx = newResIndx;
363  mass = newMass;
364  }
365 } ;
366 
367 static char *FORMAT(BigReal X)
368 {
369  static char tmp_string[25];
370  const double maxnum = 9999999999.9999;
371  if ( X > maxnum ) X = maxnum;
372  if ( X < -maxnum ) X = -maxnum;
373  sprintf(tmp_string," %14.4f",X);
374  return tmp_string;
375 }
376 
377 static char *FORMAT(const char *X)
378 {
379  static char tmp_string[25];
380  sprintf(tmp_string," %14s",X);
381  return tmp_string;
382 }
383 
384 static char *QMETITLE(int X)
385 {
386  static char tmp_string[21];
387  sprintf(tmp_string,"QMENERGY: %7d",X);
388  return tmp_string;
389 }
390 
391 
392 typedef std::pair<int, LSSDataStr> atmLSSData;
393 typedef std::map<int, LSSDataStr> LSSDataMap;
394 
395 typedef std::pair<int, Mass> refLSSData;
396 typedef std::map<int, Mass> LSSRefMap;
397 
398 typedef std::vector<ComputeQMPntChrg> QMPCVec ;
399 typedef std::vector<ComputeQMAtom> QMAtmVec ;
400 
401 class ComputeQMMgr : public CBase_ComputeQMMgr {
402 public:
403  ComputeQMMgr();
404  ~ComputeQMMgr();
405 
406  void setCompute(ComputeQM *c) { QMCompute = c; }
407 
408  // Function called on pe ZERO by all pe's.
409  // Receives partial QM coordinates form the QM atoms in each pe and aggregates
410  // them all in a single structure.
411  // This function initializes some variables and allcoates memory.
412  void recvPartQM(QMCoordMsg*) ;
413 
414  void recvFullQM(QMCoordMsg*) ;
415 
416  void recvPntChrg(QMPntChrgMsg*) ;
417 
418  void calcMOPAC(QMGrpCalcMsg *) ;
419  void calcORCA(QMGrpCalcMsg *) ;
420  void calcUSR(QMGrpCalcMsg *) ;
421 
422  void storeQMRes(QMGrpResMsg *) ;
423 
424  void procQMRes();
425 
426  void recvForce(QMForceMsg*) ;
427 
428  SortedArray<LSSSubsDat> &get_subsArray() { return subsArray; } ;
429 private:
430  ComputeQMMgr_SDAG_CODE
431 
432  CProxy_ComputeQMMgr QMProxy;
433  ComputeQM *QMCompute;
434 
435  int numSources;
436  int numArrivedQMMsg, numArrivedPntChrgMsg, numRecQMRes;
437  QMCoordMsg **qmCoordMsgs;
438  QMPntChrgMsg **pntChrgCoordMsgs;
439 
440  std::vector<int> qmPEs ;
441 
442  ComputeQMAtom *qmCoord;
443  QMPCVec pntChrgMsgVec;
444  QMForce *force;
445 
446  int numQMGrps;
447 
448  int numAtoms;
449 
450  int qmAtmIter;
451 
452  int numQMAtms;
453 
454  Bool noPC ;
455  int meNumMMIndx;
456 
457  int qmPCFreq;
458  // In case we are skiping steps between point charge re-assignment (i.e. qmPCFreq > 0),
459  // we keep a list in rank zero whihch is replaced every X steps, and
460  // then re-sent to all PEs every "X+1" steps. This is done because atoms
461  // may move between pathes and PEs between PC reassignments, if they are
462  // sufficiently long or unlucky.
463  int *pcIDList ;
464  int pcIDListSize;
465  Bool resendPCList;
467 
468  int replaceForces;
469  int bondValType;
470  SimParameters * simParams;
471  Molecule *molPtr;
472  // ID for each QM group.
473  Real *grpID;
474  Real *grpChrg, *grpMult;
475 
476  Bool qmLSSOn ;
477  int qmLSSFreq;
478  int qmLSSResSize;
479  int *qmLSSSize;
480  int *qmLSSIdxs;
481  Mass *qmLSSMass;
482  int lssTotRes;
483  // This will hold one map per QM group. Each will map an atom ID with
484  // the local residue index for the solvent molecule it belongs to,
485  // and its atomic mass (in case COM calculation is being performed).
486  ResizeArray<LSSDataMap> grpIDResNum ;
487  int *qmLSSRefIDs;
488  int *qmLSSRefSize;
489  // This will hold one map per QM group. Each will map an atom ID with
490  // its atomic mass. This is used to calculate the reference COM to which
491  // solvent molecules will be compared for distance calculation.
492  ResizeArray<LSSRefMap> lssGrpRefMass ;
493  // Current positions of each LSS RESIDUE
494  Position *lssPos;
495  Mass lssResMass;
496  SortedArray<LSSSubsDat> subsArray;
497 
498  BigReal ewaldcof, pi_ewaldcof;
499 
500  // Accumulates the enery of different QM groups. This is the energy passed to NAMD.
501  BigReal totalEnergy;
502  BigReal totVirial[3][3];
503 
504  int dcdOutFile, dcdPosOutFile;
505  Real *outputData ;
506  int timeStep ;
507 
508  void procBonds(int numBonds,
509  const int *const qmGrpBondsGrp,
510  const int *const qmMMBondedIndxGrp,
511  const int *const *const chargeTarget,
512  const int *const numTargs,
513  const QMPCVec grpPntChrgVec,
514  QMPCVec &grpAppldChrgVec) ;
515 
516  void pntChrgSwitching(QMGrpCalcMsg* msg, QMAtomData *pcPpme) ;
517 
518  void lssPrepare() ;
519  void lssUpdate(int grpIter, QMAtmVec &grpQMAtmVec, QMPCVec &grpPntChrgVec);
520 
521  void calcCSMD(int grpIndx,int numQMAtoms, const QMAtomData *atmP, Force *resForce) ;
522  Bool cSMDon;
523 
524  int cSMDnumInstances, cSMDInitGuides;
525  // Index of Conditional SMD guides per QM group
526  int const * const * cSMDindex;
527  // Num instances per QM group
528  int const * cSMDindxLen;
529  // Positions of Conditional SMD guides
530  Position* cSMDguides;
531  // Atom indices for Origin and Target of cSMD
532  int const * const * cSMDpairs;
533  // Spring constants for cSMD
534  Real const * cSMDKs;
535  // Speed of movement of guide particles for cSMD.
536  Real const * cSMDVels;
537  // Distance cutoff for guide particles for cSMD.
538  Real const * cSMDcoffs;
539  // Vector where we store all calculated cSMD forces
540  Force * cSMDForces ;
541 
542  #ifdef DEBUG_QM
543  void Write_PDB(std::string Filename, const QMGrpCalcMsg *dataMsg);
544  void Write_PDB(std::string Filename, const QMCoordMsg *dataMsg);
545  #endif
546 };
547 
549  QMProxy(thisgroup), QMCompute(0), numSources(0), numArrivedQMMsg(0),
550  numArrivedPntChrgMsg(0), numRecQMRes(0), qmCoordMsgs(0), pntChrgCoordMsgs(0),
551  qmCoord(0), force(0), numAtoms(0), dcdOutFile(0), outputData(0), pcIDListSize(0) {
552 
553  CkpvAccess(BOCclass_group).computeQMMgr = thisgroup;
554 
555 }
556 
558 
559  if (qmCoordMsgs != 0) {
560  for ( int i=0; i<numSources; ++i ) {
561  if (qmCoordMsgs[i] != 0)
562  delete qmCoordMsgs[i];
563  }
564  delete [] qmCoordMsgs;
565  }
566 
567  if (pntChrgCoordMsgs != 0) {
568  for ( int i=0; i<numSources; ++i ) {
569  if (pntChrgCoordMsgs[i] != 0)
570  delete pntChrgCoordMsgs[i];
571  }
572  delete [] pntChrgCoordMsgs;
573  }
574 
575  if (qmCoord != 0)
576  delete [] qmCoord;
577 
578  if (force != 0)
579  delete [] force;
580 
581 
582  if (dcdOutFile != 0)
583  close_dcd_write(dcdOutFile);
584 
585  if (dcdPosOutFile != 0)
586  close_dcd_write(dcdPosOutFile);
587 
588  if (outputData != 0)
589  delete [] outputData;
590 
591  if (lssPos != NULL)
592  delete [] lssPos;
593 }
594 
596  return mgr->get_subsArray();
597 } ;
598 
600  ComputeHomePatches(c), oldForces(0)
601 {
602  CProxy_ComputeQMMgr::ckLocalBranch(
603  CkpvAccess(BOCclass_group).computeQMMgr)->setCompute(this);
604 
606 
607 }
608 
610 {
611  if (oldForces != 0)
612  delete [] oldForces;
613 }
614 
616 {
618 
619  // Get some pointers that will be used in the future,
620  simParams = Node::Object()->simParameters ;
621  molPtr = Node::Object()->molecule;
622  // Future comes fast...
623  qmAtomGroup = molPtr->get_qmAtomGroup() ;
624 
625  noPC = molPtr->get_noPC();
626  if (noPC) {
627  meNumMMIndx = molPtr->get_qmMeNumBonds();
628  meMMindx = molPtr->get_qmMeMMindx();
629  meQMGrp =molPtr->get_qmMeQMGrp();
630 
631  for (int i=0; i<meNumMMIndx; i++)
632  meQMBonds.add(meMMQMGrp(meMMindx[i], meQMGrp[i]));
633  }
634 
635  numQMAtms = molPtr->get_numQMAtoms();
636  qmAtmChrg = molPtr->get_qmAtmChrg() ;
637  qmAtmIndx = molPtr->get_qmAtmIndx() ;
638 
639  numQMGrps = molPtr->get_qmNumGrps();
640  qmGrpIDArray = molPtr->get_qmGrpID() ;
641 
642  cutoff = simParams->cutoff;
643 
644  customPC = simParams->qmCustomPCSel;
645  if (customPC) {
646 
647  customPCLists.resize(numQMGrps);
648 
649  int totCustPCs = molPtr->get_qmTotCustPCs();
650  const int *custPCSizes = molPtr->get_qmCustPCSizes();
651  const int *customPCIdxs = molPtr->get_qmCustomPCIdxs();
652 
653  int minI = 0, maxI = 0, grpIter = 0;
654 
655  while (grpIter < numQMGrps) {
656 
657  maxI += custPCSizes[grpIter];
658 
659  for (size_t i=minI; i < maxI; i++) {
660 
661  customPCLists[grpIter].add( customPCIdxs[i] ) ;
662  }
663 
664  // Minimum index gets the size of the current group
665  minI += custPCSizes[grpIter];
666  // Group iterator gets incremented
667  grpIter++;
668 
669  }
670  }
671 
672 }
673 
674 
676 {
677  CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
678 
680 
681  int timeStep;
682 
683  #ifdef DEBUG_QM
684  DebugM(4,"----> Initiating QM work on rank " << CkMyPe() <<
685  " with " << patchList.size() << " patches." << std::endl );
686  #endif
687 
688  patchData.clear();
689 
690  // Determines hou many QM atoms are in this node.
691  int numLocalQMAtoms = 0, localMM1atoms = 0;
692  for (ap = ap.begin(); ap != ap.end(); ap++) {
693  CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
694  int localNumAtoms = (*ap).p->getNumAtoms();
695 
696  patchData.push_back(patchDataStrc((*ap).positionBox,
697  (*ap).positionBox->open(),
698  (*ap).p) );
699 
700  for(int i=0; i<localNumAtoms; ++i) {
701  if ( qmAtomGroup[xExt[i].id] > 0 ) {
702  numLocalQMAtoms++;
703  }
704  else if (meNumMMIndx > 0) {
705  // If we have a mechanical embedding simulation with no
706  // point charges, but with QM-MM bonds, we look for the MM1 atoms
707  // here and pass them directly. No need for an extra round of
708  // comms just to get atoms we already know we need.
709 
710  auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
711 
712  if (retIt != NULL) {
713  localMM1atoms++;
714  }
715  }
716  }
717 
718  timeStep = (*ap).p->flags.step ;
719  }
720 
721  DebugM(4, "Node " << CkMyPe() << " has " << numLocalQMAtoms
722  << " QM atoms." << std::endl) ;
723  #ifdef DEBUG_QM
724  if (localMM1atoms > 0)
725  DebugM(4, "Node " << CkMyPe() << " has " << localMM1atoms
726  << " MM1 atoms." << std::endl) ;
727  #endif
728  // Send QM atoms
729 
730  // This pointer will be freed in rank zero.
731  QMCoordMsg *msg = new (numLocalQMAtoms + localMM1atoms,0, 0) QMCoordMsg;
732  msg->sourceNode = CkMyPe();
733  msg->numAtoms = numLocalQMAtoms + localMM1atoms;
734  msg->timestep = timeStep;
735  ComputeQMAtom *data_ptr = msg->coord;
736 
737  // Build a message with coordinates, charge, atom ID and QM group
738  // for all QM atoms in the node, then send it to node zero.
739  int homeCounter = 0;
740  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
741 
742  CompAtom *x = (*pdIt).compAtomP;
743  const CompAtomExt *xExt = (*pdIt).homePatchP->getCompAtomExtInfo();
744  int localNumAtoms = (*pdIt).homePatchP->getNumAtoms();
745  const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
746 
747  for(int i=0; i<localNumAtoms; ++i) {
748 
749  if ( qmAtomGroup[xExt[i].id] > 0 ) {
750 
751  Real charge = 0;
752 
753  for (int qi=0; qi<numQMAtms; qi++) {
754  if (qmAtmIndx[qi] == xExt[i].id) {
755  charge = qmAtmChrg[qi];
756  break;
757  }
758  }
759 
760  data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position,
761  fullAtms[i].transform) ;
762 // data_ptr->charge = x[i].charge;
763  data_ptr->charge = charge;
764  data_ptr->id = xExt[i].id;
765  data_ptr->qmGrpID = qmAtomGroup[xExt[i].id] ;
766  data_ptr->homeIndx = homeCounter;
767  data_ptr->vdwType = fullAtms[i].vdwType;
768 
769  ++data_ptr;
770  }
771  else if (meNumMMIndx > 0) {
772  // If we have a mechanical embedding simulation with no
773  // point charges, but with QM-MM bonds, we look for the MM1 atoms
774  // connected here and pass them directly.
775 
776  auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
777 
778  if (retIt != NULL) {
779 
780  DebugM(3,"Found atom " << retIt->mmIndx << "," << retIt->qmGrp << "\n" );
781 
782  data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position,
783  fullAtms[i].transform) ;
784  // We force the charge to be zero since we would only get
785  // here if mechanical embeding was selected.
786  data_ptr->charge = 0;
787  data_ptr->id = xExt[i].id;
788  data_ptr->qmGrpID = retIt->qmGrp ;
789  data_ptr->homeIndx = homeCounter;
790 
791  // We re-use this varaible to indicate this is an MM1 atom.
792  data_ptr->vdwType = -1;
793 
794  ++data_ptr;
795 
796  }
797 
798  }
799 
800  homeCounter++;
801 
802  }
803 
804  }
805 
806  if (noPC) {
807  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
808  CompAtom *x = (*pdIt).compAtomP;
809  (*pdIt).posBoxP->close(&x);
810  }
811  }
812 
813  QMProxy[0].recvPartQM(msg);
814 
815 }
816 
817 
819 {
820  // In the first (ever) step of the simulation, we allocate arrays that
821  // will be used throughout the simulation, and get some important numbers.
822  if ( ! numSources ) {
823  DebugM(4,"Initializing ComputeQMMgr variables." << std::endl);
824  numSources = (PatchMap::Object())->numNodesWithPatches();
825 
826  DebugM(4,"There are " << numSources << " nodes with patches." << std::endl);
827  qmCoordMsgs = new QMCoordMsg*[numSources];
828  for ( int i=0; i<numSources; ++i ) {
829  qmCoordMsgs[i] = 0;
830  }
831 
832  // Prepares the allocation for the recvPntChrg function.
833 
834  DebugM(4,"Getting data from molecule and simParameters." << std::endl);
835 
836  molPtr = Node::Object()->molecule;
837  simParams = Node::Object()->simParameters;
838 
839  numAtoms = molPtr->numAtoms;
840  force = new QMForce[numAtoms];
841 
842  numQMAtms = molPtr->get_numQMAtoms();
843  qmAtmIter = 0;
844 
845  noPC = simParams->qmNoPC ;
846  meNumMMIndx = molPtr->get_qmMeNumBonds();
847  if (noPC && meNumMMIndx == 0) {
848  pntChrgCoordMsgs = NULL;
849  }
850  else {
851  pntChrgCoordMsgs = new QMPntChrgMsg*[numSources];
852  for ( int i=0; i<numSources; ++i ) {
853  pntChrgCoordMsgs[i] = 0;
854  }
855  }
856 
857  qmPCFreq = molPtr->get_qmPCFreq();
858  resendPCList = false;
859 
860  grpID = molPtr->get_qmGrpID() ;
861  bondValType = simParams->qmBondValType;
862 
863  numQMGrps = molPtr->get_qmNumGrps();
864 
865  grpChrg = molPtr->get_qmGrpChrg() ;
866 
867  grpMult = molPtr->get_qmGrpMult() ;
868 
869  qmLSSOn = simParams->qmLSSOn ;
870  if (qmLSSOn) {
871  qmLSSFreq = molPtr->get_qmLSSFreq() ;
872  qmLSSSize = molPtr->get_qmLSSSize() ;
873  qmLSSIdxs = molPtr->get_qmLSSIdxs() ;
874  qmLSSMass = molPtr->get_qmLSSMass() ;
875  qmLSSResSize = molPtr->get_qmLSSResSize() ;
876  qmLSSRefIDs = molPtr->get_qmLSSRefIDs() ;
877  qmLSSRefSize = molPtr->get_qmLSSRefSize() ;
878 
879  lssPrepare();
880  }
881 
882  numArrivedQMMsg = 0 ;
883  numArrivedPntChrgMsg = 0 ;
884 
885  qmCoord = new ComputeQMAtom[numQMAtms];
886 
887  replaceForces = 0;
888  if (molPtr->get_qmReplaceAll()) {
889  replaceForces = 1;
890  }
891 
892  cSMDon = molPtr->get_qmcSMD() ;
893  if (cSMDon) {
894 
895  // We have to initialize the guide particles during the first step.
896  cSMDInitGuides = 0;
897 
898  cSMDnumInstances = molPtr->get_cSMDnumInst();
899  cSMDindex = molPtr->get_cSMDindex();
900  cSMDindxLen = molPtr->get_cSMDindxLen();
901  cSMDpairs = molPtr->get_cSMDpairs();
902  cSMDKs = molPtr->get_cSMDKs();
903  cSMDVels = molPtr->get_cSMDVels();
904  cSMDcoffs = molPtr->get_cSMDcoffs();
905 
906  cSMDguides = new Position[cSMDnumInstances];
907  cSMDForces = new Force[cSMDnumInstances];
908  }
909 
910 
911  DebugM(4,"Initializing DCD file for charge information." << std::endl);
912 
913  // Initializes output DCD file for charge information.
914  if (simParams->qmOutFreq > 0) {
915  std::string dcdOutFileStr;
916  dcdOutFileStr.clear();
917  dcdOutFileStr.append(simParams->outputFilename) ;
918  dcdOutFileStr.append(".qdcd") ;
919  dcdOutFile = open_dcd_write(dcdOutFileStr.c_str()) ;
920 
921  if (dcdOutFile == DCD_FILEEXISTS) {
922  iout << iERROR << "DCD file " << dcdOutFile << " already exists!!\n" << endi;
923  NAMD_err("Could not write QM charge DCD file.");
924  } else if (dcdOutFile < 0) {
925  iout << iERROR << "Couldn't open DCD file " << dcdOutFile << ".\n" << endi;
926  NAMD_err("Could not write QM charge DCD file.");
927  } else if (! dcdOutFile) {
928  NAMD_bug("ComputeQMMgr::recvPartQM open_dcd_write returned fileid of zero");
929  }
930 
931  timeStep = simParams->firstTimestep;
932 
933  int NSAVC, NFILE, NPRIV, NSTEP;
934  NSAVC = simParams->qmOutFreq;
935  NPRIV = timeStep;
936  NSTEP = NPRIV - NSAVC;
937  NFILE = 0;
938 
939  // Write out the header
940  int ret_code = write_dcdheader(dcdOutFile, dcdOutFileStr.c_str(),
941  numQMAtms, NFILE, NPRIV, NSAVC, NSTEP,
942  simParams->dt/TIMEFACTOR, 0);
943 
944  if (ret_code<0) {
945  NAMD_err("Writing of DCD header failed!!");
946  }
947 
948  // The DCD write function takes 3 independent arrays for X, Y and Z
949  // coordinates, but we allocate one and send in the pieces.
950  outputData = new Real[3*numQMAtms];
951  }
952 
953  DebugM(4,"Initializing DCD file for position information." << std::endl);
954  // Initializes output DCD file for position information.
955  if (simParams->qmPosOutFreq > 0) {
956  std::string dcdPosOutFileStr;
957  dcdPosOutFileStr.clear();
958  dcdPosOutFileStr.append(simParams->outputFilename) ;
959  dcdPosOutFileStr.append(".QMonly.dcd") ;
960  dcdPosOutFile = open_dcd_write(dcdPosOutFileStr.c_str()) ;
961 
962  if (dcdPosOutFile == DCD_FILEEXISTS) {
963  iout << iERROR << "DCD file " << dcdPosOutFile << " already exists!!\n" << endi;
964  NAMD_err("Could not write QM charge DCD file.");
965  } else if (dcdPosOutFile < 0) {
966  iout << iERROR << "Couldn't open DCD file " << dcdPosOutFile << ".\n" << endi;
967  NAMD_err("Could not write QM charge DCD file.");
968  } else if (! dcdPosOutFile) {
969  NAMD_bug("ComputeQMMgr::recvPartQM open_dcd_write returned fileid of zero");
970  }
971 
972  timeStep = simParams->firstTimestep;
973 
974  int NSAVC, NFILE, NPRIV, NSTEP;
975  NSAVC = simParams->qmOutFreq;
976  NPRIV = timeStep;
977  NSTEP = NPRIV - NSAVC;
978  NFILE = 0;
979 
980  // Write out the header
981  int ret_code = write_dcdheader(dcdPosOutFile, dcdPosOutFileStr.c_str(),
982  numQMAtms, NFILE, NPRIV, NSAVC, NSTEP,
983  simParams->dt/TIMEFACTOR, 0);
984 
985  if (ret_code<0) {
986  NAMD_err("Writing of DCD header failed!!");
987  }
988 
989  // The DCD write function takes 3 independent arrays for X, Y and Z
990  // coordinates, but we allocate one and send in the pieces.
991  outputData = new Real[3*numQMAtms];
992  }
993 
994  // Prepares list of PEs which will run the QM software
995  int simsPerNode = simParams->qmSimsPerNode ;
996  int zeroNodeSize = CmiNumPesOnPhysicalNode(0);
997 
998  // Check if the node has enought PEs to run the requested number of simulations.
999  if ( zeroNodeSize < simsPerNode ) {
1000  iout << iERROR << "There are " << zeroNodeSize << " cores pernode, but "
1001  << simsPerNode << " QM simulations per node were requested.\n" << endi ;
1002  NAMD_die("Error preparing QM simulations.");
1003  }
1004 
1005  DebugM(4,"Preparing PE list for QM execution.\n");
1006  qmPEs.clear(); // Making sure its empty.
1007 
1008  int numNodes = CmiNumPhysicalNodes();
1009  int numPlacedQMGrps = 0;
1010  int placedOnNode = 0;
1011  int nodeIt = 0 ;
1012 
1013  // The default is to only run on rank zero.
1014  if ( simsPerNode <= 0 ) {
1015  qmPEs.push_back(0);
1016  numPlacedQMGrps = 1;
1017  }
1018 
1019  while ( (numPlacedQMGrps < numQMGrps) && (simsPerNode > 0) ) {
1020 
1021  // If we searched all nodes, break the loop.
1022  if (nodeIt == numNodes) {
1023  break;
1024  }
1025 
1026  int *pelist;
1027  int nodeSize;
1028  CmiGetPesOnPhysicalNode(CmiPhysicalNodeID(nodeIt), &pelist, &nodeSize);
1029 
1030  DebugM(4,"Checking node " << nodeIt +1 << " out of " << numNodes
1031  << " (" << nodeSize << " PEs: " << pelist[0] << " to "
1032  << pelist[nodeSize-1] << ")." << std::endl );
1033 
1034  for ( int i=0; i<nodeSize; ++i ) {
1035 
1036  // Check if the PE has patches. We only run on PEs with patches!
1037  if ( (PatchMap::Object())->numPatchesOnNode(pelist[i]) == 0 ) {
1038  DebugM(1,"PE " << pelist[i] << " has no patches! We'll skip it."
1039  << std::endl);
1040  continue ;
1041  }
1042 
1043  // Add the target PE on the target node to the list
1044  // of PEs which will carry QM simulations.
1045  qmPEs.push_back(pelist[i]);
1046 
1047  DebugM(1,"Added PE to QM execution list: " << pelist[i] << "\n");
1048 
1049  numPlacedQMGrps++;
1050  placedOnNode++;
1051 
1052  if (placedOnNode == simsPerNode) {
1053  DebugM(1,"Placed enought simulations on this node.\n");
1054  break;
1055  }
1056 
1057 
1058  }
1059 
1060  nodeIt++ ;
1061  placedOnNode = 0;
1062  }
1063 
1064  if ( numPlacedQMGrps < numQMGrps ) {
1065  iout << iWARN << "Could not compute all QM groups in parallel.\n" << endi ;
1066  }
1067 
1068  iout << iINFO << "List of ranks running QM simulations: " << qmPEs[0] ;
1069  for (int i=1; i < qmPEs.size(); i++) {
1070  iout << ", " << qmPEs[i] ;
1071  }
1072  iout << ".\n" << endi;
1073 
1074  }
1075 
1076  DebugM(1,"Receiving from rank " << msg->sourceNode
1077  << " a total of " << msg->numAtoms << " QM atoms." << std::endl);
1078 
1079  // In case we are NOT using point charges but there are QM-MM bonds,
1080  // test each QM message for MM1 atoms.
1081  if (meNumMMIndx > 0) {
1082 
1084  ComputeQMAtom *data_ptr = msg->coord;
1085 
1086  for (int i=0; i<msg->numAtoms; i++) {
1087  if (data_ptr[i].vdwType < 0) {
1088  tmpAtm.add(data_ptr[i]) ;
1089  }
1090  }
1091 
1092  QMPntChrgMsg *pntChrgMsg = new (tmpAtm.size(), 0) QMPntChrgMsg;
1093  pntChrgMsg->sourceNode = msg->sourceNode ;
1094  pntChrgMsg->numAtoms = tmpAtm.size() ;
1095  ComputeQMPntChrg* newPCData = pntChrgMsg->coord ;
1096 
1097  QMCoordMsg *newMsg = msg;
1098 
1099  if (tmpAtm.size() > 0) {
1100 
1101  newMsg = new (msg->numAtoms - tmpAtm.size(),0, 0) QMCoordMsg;
1102  newMsg->sourceNode = msg->sourceNode ;
1103  newMsg->numAtoms = msg->numAtoms - tmpAtm.size() ;
1104  newMsg->timestep = msg->timestep ;
1105  ComputeQMAtom *newMsgData = newMsg->coord;
1106 
1107  for (int i=0; i<msg->numAtoms; i++) {
1108  if (data_ptr[i].vdwType < 0) {
1109  newPCData->position = data_ptr[i].position ;
1110  newPCData->charge = data_ptr[i].charge ;
1111  newPCData->id = data_ptr[i].id ;
1112  newPCData->qmGrpID = data_ptr[i].qmGrpID ;
1113  newPCData->homeIndx = data_ptr[i].homeIndx ;
1114  newPCData->dist = 0 ;
1115  newPCData->mass = 0 ;
1116  newPCData->vdwType = 0 ;
1117  newPCData++;
1118  }
1119  else {
1120  *newMsgData = data_ptr[i] ;
1121  newMsgData++;
1122  }
1123  }
1124 
1125  delete msg;
1126 
1127  }
1128 
1129  qmCoordMsgs[numArrivedQMMsg] = newMsg;
1130  ++numArrivedQMMsg;
1131 
1132  pntChrgCoordMsgs[numArrivedPntChrgMsg] = pntChrgMsg;
1133  ++numArrivedPntChrgMsg;
1134  }
1135  else {
1136  qmCoordMsgs[numArrivedQMMsg] = msg;
1137  ++numArrivedQMMsg;
1138  }
1139 
1140  if ( numArrivedQMMsg < numSources )
1141  return;
1142 
1143  // Now that all messages arrived, get all QM positions.
1144  for (int msgIt=0; msgIt < numArrivedQMMsg; msgIt++){
1145 
1146  DebugM(1, "Getting positions for " << qmCoordMsgs[msgIt]->numAtoms
1147  << " QM atoms in this message." << std::endl);
1148 
1149  for ( int i=0; i < qmCoordMsgs[msgIt]->numAtoms; ++i ) {
1150  qmCoord[qmAtmIter] = qmCoordMsgs[msgIt]->coord[i];
1151  qmAtmIter++;
1152  }
1153  }
1154 
1155  if (qmAtmIter != numQMAtms) {
1156  iout << iERROR << "The number of QM atoms received (" << qmAtmIter
1157  << ") is different than expected: " << numQMAtms << "\n" << endi;
1158  NAMD_err("Problems broadcasting QM atoms.");
1159  }
1160 
1161  // Resets the counter for the next step.
1162  numArrivedQMMsg = 0;
1163  qmAtmIter = 0;
1164 
1165  timeStep = qmCoordMsgs[0]->timestep;
1166 
1167  // Makes sure there is no trash or old info in the force array.
1168  // This is where we accumulate forces from the QM software and our own
1169  // Coulomb forces. It will have info on QM atoms and point charges only.
1170  for (int i=0; i<numAtoms; ++i ) {
1171  force[i].force = 0; // this assigns 0 (ZERO) to all componenets of the vector.
1172  force[i].replace = 0; // By default, no atom has all forces replaced.
1173  force[i].homeIndx = -1; // To prevent errors from sliping.
1174  force[i].charge = 0;
1175  force[i].id = i; // Initializes the variable since not all forces are sent to all ranks.
1176  }
1177 
1178 
1179  for (int i=0; i<numQMAtms; i++) {
1180  // Each force receives the home index of its atom with respect to the
1181  // local set of atoms in each node.
1182  if (force[qmCoord[i].id].homeIndx != -1
1183  && force[qmCoord[i].id].homeIndx != qmCoord[i].homeIndx
1184  ) {
1185  iout << iERROR << "Overloading QM atom "
1186  << qmCoord[i].id << "; home index: "
1187  << force[qmCoord[i].id].homeIndx << " with " << qmCoord[i].homeIndx
1188  << "\n" << endi ;
1189  NAMD_die("Error preparing QM simulations.");
1190  }
1191 
1192  force[qmCoord[i].id].homeIndx = qmCoord[i].homeIndx;
1193  // Each force on QM atoms has an indicator so NAMD knows if all
1194  // NAMD forces should be erased and completely overwritten by the
1195  // external QM forces.
1196  force[qmCoord[i].id].replace = replaceForces;
1197  }
1198 
1199  if (noPC) {
1200  // this pointer should be freed in rank zero, after receiving it.
1201  QMPntChrgMsg *pntChrgMsg = new (0, 0) QMPntChrgMsg;
1202  pntChrgMsg->sourceNode = CkMyPe();
1203  pntChrgMsg->numAtoms = 0;
1204 
1205  QMProxy[0].recvPntChrg(pntChrgMsg);
1206  }
1207  else {
1208  // The default mode is to update the poitn charge selection at every step.
1209  int pcSelMode = PCMODEUPDATESEL;
1210 
1211  int msgPCListSize = 0;
1212  // We check wether point charges are to be selected with a stride.
1213  if ( qmPCFreq > 0 ) {
1214  if (timeStep % qmPCFreq == 0) {
1215  // Marks that the PC list determined in this step will
1216  // be broadcast on the *next* step, and kept for the following
1217  // qmPCFreq-1 steps.
1218  resendPCList = true;
1219 
1220  // Clears the list since we don't know how many charges
1221  // will be selected this time.
1222  delete [] pcIDList;
1223  }
1224  else {
1225  // If the PC selection is not to be updated in this step,
1226  // inform the nodes that the previous list (or the updated
1227  // list being sent in this message) is to be used and only
1228  // updated positions will be returned.
1229  pcSelMode = PCMODEUPDATEPOS;
1230 
1231  // If this is the first step after a PC re-selection, all
1232  // ranks receive the total list, since atoms may move between
1233  // PEs in between PC re-assignments (if they are far enought apart
1234  // or if you are unlucky)
1235  if (resendPCList) {
1236  msgPCListSize = pcIDListSize;
1237  resendPCList = false;
1238  }
1239  }
1240  }
1241 
1242  // In case we are using custom selection of point charges, indicate the mode.
1243  if (simParams->qmCustomPCSel)
1244  pcSelMode = PCMODECUSTOMSEL;
1245 
1246  DebugM(1,"Broadcasting current positions of QM atoms.\n");
1247  for ( int j=0; j < numSources; ++j ) {
1248  // This pointer will be freed in the receiving rank.
1249  QMCoordMsg *qmFullMsg = new (numQMAtms, msgPCListSize, 0) QMCoordMsg;
1250  qmFullMsg->sourceNode = CkMyPe();
1251  qmFullMsg->numAtoms = numQMAtms;
1252  qmFullMsg->pcSelMode = pcSelMode;
1253  qmFullMsg->numPCIndxs = msgPCListSize;
1254  ComputeQMAtom *data_ptr = qmFullMsg->coord;
1255  int *msgPCListP = qmFullMsg->pcIndxs;
1256 
1257  for (int i=0; i<numQMAtms; i++) {
1258  data_ptr->position = qmCoord[i].position;
1259  data_ptr->charge = qmCoord[i].charge;
1260  data_ptr->id = qmCoord[i].id;
1261  data_ptr->qmGrpID = qmCoord[i].qmGrpID;
1262  data_ptr->homeIndx = -1; // So there is no mistake that there is no info here.
1263  data_ptr++;
1264  }
1265 
1266  for (int i=0; i<msgPCListSize; i++) {
1267  msgPCListP[i] = pcIDList[i] ;
1268  }
1269 
1270  #ifdef DEBUG_QM
1271  if (j == 0)
1272  Write_PDB("qmMsg.pdb", qmFullMsg) ;
1273  #endif
1274 
1275  // The messages are deleted later, we will need them.
1276  QMProxy[qmCoordMsgs[j]->sourceNode].recvFullQM(qmFullMsg);
1277  }
1278  }
1279 }
1280 
1282 
1283  if (subsArray.size() > 0)
1284  subsArray.clear();
1285 
1286  QMCompute->processFullQM(qmFullMsg);
1287 }
1288 
1289 typedef std::map<Real, std::pair<Position,BigReal> > GrpDistMap ;
1290 typedef std::pair<Position,BigReal> PosDistPair ;
1291 
1293 
1295 
1296  const Lattice & lattice = patchList[0].p->lattice;
1297 
1298  // Dynamically accumulates point charges as they are found.
1299  // The same MM atom may be added to this vector more than once if it is
1300  // within the cutoff region of two or more QM regions.
1301  ResizeArray<ComputeQMPntChrg> compPCVec ;
1302 
1303  // This will keep the set of QM groups with which each
1304  // point charge should be used. It is re-used for each atom in this
1305  // patch so we can controll the creation of the compPCVec vector.
1306  // The first item in the pair is the QM Group ID, the second is a pair
1307  // with the point charge position (after PBC wraping) and the shortest
1308  // distance between the point charge and an atom of this QM group.
1309  GrpDistMap chrgGrpMap ;
1310 
1311  DebugM(4,"Rank " << CkMyPe() << " receiving from rank " << qmFullMsg->sourceNode
1312  << " a total of " << qmFullMsg->numAtoms << " QM atoms and "
1313  << qmFullMsg->numPCIndxs << " PC IDs." << std::endl);
1314 
1315  // If this is the firts step after a step of PC re-selection,
1316  // we store all PC IDs that all PEs found, in case some atoms end up
1317  // here untill the next re-selection step.
1318  if (qmFullMsg->numPCIndxs) {
1319 
1320  pcIDSortList.clear();
1321 
1322  int *pcIndx = qmFullMsg->pcIndxs;
1323  for (int i=0; i< qmFullMsg->numPCIndxs;i++) {
1324  pcIDSortList.load(pcIndx[i]);
1325  }
1326 
1327  pcIDSortList.sort();
1328  }
1329 
1330  int totalNodeAtoms = 0;
1331  int atomIter = 0;
1332  int uniqueCounter = 0;
1333 
1334  const ComputeQMAtom *qmDataPtn = qmFullMsg->coord;
1335 
1336  switch ( qmFullMsg->pcSelMode ) {
1337 
1338  case PCMODEUPDATESEL:
1339  {
1340 
1341  DebugM(4,"Updating PC selection.\n")
1342 
1343  // Loops over all atoms in this node and checks if any MM atom is within
1344  // the cutof radius from a QM atom
1345  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
1346 
1347  CompAtom *x = (*pdIt).compAtomP;
1348  const CompAtomExt *xExt = (*pdIt).homePatchP->getCompAtomExtInfo();
1349  int localNumAtoms = (*pdIt).homePatchP->getNumAtoms();
1350  const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
1351 
1352  totalNodeAtoms += localNumAtoms;
1353 
1354  // Iterates over the local atoms in a PE, all patches.
1355  for(int i=0; i<localNumAtoms; ++i) {
1356 
1357  // A new subset for each atom in this node.
1358  chrgGrpMap.clear();
1359 
1360  Real pcGrpID = qmAtomGroup[xExt[i].id];
1361 
1362  Real charge = x[i].charge;
1363 
1364  // If the atom is a QM atom, there will be no charge info in the
1365  // atom box. Therefore, we take the charge in the previous step
1366  // in case this atom is a point charge for a neighboring QM region.
1367  if (pcGrpID > 0) {
1368  for (int qi=0; qi<numQMAtms; qi++) {
1369  if (qmAtmIndx[qi] == xExt[i].id) {
1370  charge = qmAtmChrg[qi];
1371  break;
1372  }
1373  }
1374  }
1375 
1376  qmDataPtn = qmFullMsg->coord;
1377  for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
1378 
1379  Real qmGrpID = qmDataPtn->qmGrpID;
1380 
1381  // We check If a QM atom and the node atom "i"
1382  // belong to the same QM group. If not, atom "i" is either an MM
1383  // atom or a QM atom from another group, which will be seen
1384  // as an MM point charge.
1385  // If they belong to the same group, just skip the distance
1386  // calculation and move on to the next QM atom.
1387  // This loop needs to go over all QM atoms since a point charge
1388  // may me be close to two different QM groups, in which case it
1389  // will be sent to both.
1390  if (qmGrpID == pcGrpID) {
1391  continue;
1392  }
1393 
1394  // calculate distance between QM atom and Poitn Charge
1395  Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
1396 
1397  BigReal dist = r12.length(); // Distance between atoms
1398 
1399  if ( dist < cutoff ){
1400 
1401  // Since the QM region may have been wrapped around the periodic box, we
1402  // reposition point charges with respect to unwrapped coordinates of QM atoms,
1403  // finding the shortest distance between the point charge and the QM region.
1404 
1405  Position posMM = qmDataPtn->position + r12;
1406  auto ret = chrgGrpMap.find(qmGrpID) ;
1407 
1408  // If 'ret' points to end, it means the item was not found
1409  // and will be added, which means a new QM region has this
1410  // atom within its cutoff region,
1411  // which means a new point charge will be
1412  // created in the QM system in question.
1413  // 'ret' means a lot to me!
1414  if ( ret == chrgGrpMap.end()) {
1415  chrgGrpMap.insert(std::pair<Real,PosDistPair>(qmGrpID,
1416  PosDistPair(posMM,dist) ) );
1417  }
1418  else {
1419  // If the current distance is smaller than the one in
1420  // the map, we update it so that we have the smallest
1421  // distance between the point charge and any QM atom in
1422  // this group.
1423  if (dist < ret->second.second) {
1424  ret->second.first = posMM ;
1425  ret->second.second = dist ;
1426  }
1427  }
1428  }
1429  }
1430 
1431  for(auto mapIt = chrgGrpMap.begin();
1432  mapIt != chrgGrpMap.end(); mapIt++) {
1433 
1434  // We now add the final info about this point charge
1435  // to the vector, repeating it for each QM group that has it
1436  // within its cuttoff radius.
1437  compPCVec.add(
1438  ComputeQMPntChrg(mapIt->second.first, charge, xExt[i].id,
1439  mapIt->first, atomIter, mapIt->second.second,
1440  fullAtms[i].mass, fullAtms[i].vdwType)
1441  );
1442  }
1443 
1444  // Counts how many atoms are seens as point charges, by one or more
1445  // QM groups.
1446  if (chrgGrpMap.size() > 0)
1447  uniqueCounter++;
1448 
1449  atomIter++;
1450  }
1451 
1452  (*pdIt).posBoxP->close(&x);
1453  }
1454 
1455  }break; // End case PCMODEUPDATESEL
1456 
1457  case PCMODEUPDATEPOS:
1458  {
1459 
1460  DebugM(4,"Updating PC positions ONLY.\n")
1461 
1462  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
1463 
1464  CompAtom *x = (*pdIt).compAtomP;
1465  const CompAtomExt *xExt = (*pdIt).homePatchP->getCompAtomExtInfo();
1466  int localNumAtoms = (*pdIt).homePatchP->getNumAtoms();
1467  const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
1468 
1469  totalNodeAtoms += localNumAtoms;
1470 
1471  // Iterates over the local atoms in a PE, all patches.
1472  for(int i=0; i<localNumAtoms; ++i) {
1473 
1474  if (pcIDSortList.find(xExt[i].id) != NULL ) {
1475 
1476  BigReal dist = 0;
1477  Position posMM = 0; // Temporary point charge position for QM calculation.
1478 
1479  bool firstIngroupQM = true;
1480  qmDataPtn = qmFullMsg->coord;
1481  for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
1482 
1483  // Assuming we only have one QM region for now
1484  // This optio will be deprecated
1485 
1486  // Only verify distances to atoms in the QM group for which we are
1487  // gathering point charges.
1488 // if ( qmDataPtn->qmGrpID != qmGrpIDArray[grpIndx] ) continue;
1489 
1490  // calculate distance between QM atom and Poitn Charge
1491  Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
1492 
1493  // We wait to find the first QM atom in the target QM group before we
1494  // set the first guess at a minimal distance and point charge position.
1495  if ( firstIngroupQM ) {
1496  dist = r12.length();
1497 
1498  // Reposition classical point charge with respect to unwrapped QM position.
1499  posMM = qmDataPtn->position + r12;
1500 
1501  firstIngroupQM = false;
1502  }
1503 
1504  // If we find a QM atom with a shorter dstance to the point charge,
1505  // set that as the PC distance and reposition the PC
1506  if ( r12.length() < dist ) {
1507  dist = r12.length();
1508  posMM = qmDataPtn->position + r12 ;
1509  }
1510 
1511  }
1512 
1513  Real pcGrpID = qmAtomGroup[xExt[i].id];
1514  Real charge = x[i].charge;
1515 
1516  // If the atom is a QM atom, there will be no charge info in the
1517  // atom box. Therefore, we take the charge in the previous step
1518  // in case this atom is a point charge for a neighboring QM region.
1519  if (pcGrpID > 0) {
1520  for (int qi=0; qi<numQMAtms; qi++) {
1521  if (qmAtmIndx[qi] == xExt[i].id) {
1522  charge = qmAtmChrg[qi];
1523  break;
1524  }
1525 
1526  }
1527  }
1528 
1529  compPCVec.add(
1530  ComputeQMPntChrg(posMM, charge, xExt[i].id,
1531  0, atomIter, 0,
1532  fullAtms[i].mass, fullAtms[i].vdwType));
1533  }
1534 
1535  atomIter++;
1536  }
1537 
1538  (*pdIt).posBoxP->close(&x);
1539  }
1540  }break ; // End case PCMODEUPDATEPOS
1541 
1542  case PCMODECUSTOMSEL:
1543  {
1544 
1545  DebugM(4,"Updating PC positions for custom PC selection.\n")
1546 
1547  for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
1548 
1549  CompAtom *x = (*pdIt).compAtomP;
1550  const CompAtomExt *xExt = (*pdIt).homePatchP->getCompAtomExtInfo();
1551  int localNumAtoms = (*pdIt).homePatchP->getNumAtoms();
1552  const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
1553 
1554  totalNodeAtoms += localNumAtoms;
1555 
1556  // Iterates over the local atoms in a PE, all patches.
1557  for(int i=0; i<localNumAtoms; ++i) {
1558 
1559  // Checks if the atom ID belongs to a point charge list,
1560  // which indicates it is to be sent to rank zero for QM calculations.
1561  // With one list per QM group, the same point charge can be added to
1562  // different QM groups, as if it was being dynamically selected.
1563  for (int grpIndx=0; grpIndx<numQMGrps; grpIndx++) {
1564 
1565  if (customPCLists[grpIndx].find(xExt[i].id) != NULL){
1566 
1567  // Initialize the search for the smallest distance between
1568  // point charge and QM atoms. This will determine the closest position
1569  // of the point charge accounting for periodic boundary conditions.
1570 
1571  // Since the QM region may have been wrapped around the periodic box, we
1572  // reposition point charges with respect to unwrapped coordinates of QM atoms.
1573 
1574  BigReal dist = 0;
1575  Position posMM = 0; // Temporary point charge position for QM calculation.
1576 
1577  bool firstIngroupQM = true;
1578  qmDataPtn = qmFullMsg->coord;
1579  for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
1580 
1581  // Only verify distances to atoms in the QM group for which we are
1582  // gathering point charges.
1583  if ( qmDataPtn->qmGrpID != qmGrpIDArray[grpIndx] ) continue;
1584 
1585  // calculate distance between QM atom and Poitn Charge
1586  Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
1587 
1588  // We wait to find the first QM atom in the target QM group before we
1589  // set the first guess at a minimal distance and point charge position.
1590  if ( firstIngroupQM ) {
1591  dist = r12.length();
1592 
1593  // Reposition classical point charge with respect to unwrapped QM position.
1594  posMM = qmDataPtn->position + r12;
1595 
1596  firstIngroupQM = false;
1597  }
1598 
1599  // If we find a QM atom with a shorter dstance to the point charge,
1600  // set that as the PC distance and reposition the PC
1601  if ( r12.length() < dist ) {
1602  dist = r12.length();
1603  posMM = qmDataPtn->position + r12 ;
1604  }
1605 
1606  }
1607 
1608  // After determining the position o fthe point charge, get its partial charge.
1609  Real pcGrpID = qmAtomGroup[xExt[i].id];
1610  Real charge = x[i].charge;
1611 
1612  // If the atom is a QM atom, there will be no charge info in the
1613  // atom box. Therefore, we take the charge in the previous step
1614  // in case this atom is a point charge for a neighboring QM region.
1615  if (pcGrpID > 0) {
1616  for (int qi=0; qi<numQMAtms; qi++) {
1617  if (qmAtmIndx[qi] == xExt[i].id) {
1618  charge = qmAtmChrg[qi];
1619  break;
1620  }
1621  }
1622  }
1623 
1624  compPCVec.add(
1625  ComputeQMPntChrg(posMM, charge, xExt[i].id,
1626  qmGrpIDArray[grpIndx], atomIter, dist,
1627  fullAtms[i].mass, fullAtms[i].vdwType));
1628  }
1629 
1630  }
1631 
1632  atomIter++;
1633  }
1634 
1635  (*pdIt).posBoxP->close(&x);
1636  }
1637 
1638  } break;
1639  }
1640 
1641  DebugM(4,"Rank " << CkMyPe() << " found a total of " << compPCVec.size()
1642  << " point charges, out of " << totalNodeAtoms
1643  << " atoms in this node. " << uniqueCounter << " are unique." << std::endl);
1644 
1645  // Send only the MM atoms within radius
1646 
1647  // this pointer should be freed in rank zero, after receiving it.
1648  QMPntChrgMsg *pntChrgMsg = new (compPCVec.size(), 0) QMPntChrgMsg;
1649  pntChrgMsg->sourceNode = CkMyPe();
1650  pntChrgMsg->numAtoms = compPCVec.size();
1651 
1652  for (int i=0; i<compPCVec.size(); i++ ) {
1653 
1654  // Only sends the positions and charges of atoms within
1655  // cutoff of a (any) QM atom. Repeats for the number of QM groups
1656  // this charge is near to, and indicates the QM group it should
1657  // be used with.
1658  pntChrgMsg->coord[i] = compPCVec[i] ;
1659 
1660  }
1661 
1662  DebugM(4,"Rank " << pntChrgMsg->sourceNode << " sending a total of "
1663  << compPCVec.size() << " elements to rank zero." << std::endl);
1664 
1665  CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
1666  QMProxy[0].recvPntChrg(pntChrgMsg);
1667 
1668  delete qmFullMsg;
1669 }
1670 
1671 
1672 
1673 void ComputeQMMgr::procBonds(int numBonds,
1674  const int *const qmGrpBondsGrp,
1675  const int *const qmMMBondedIndxGrp,
1676  const int *const *const chargeTarget,
1677  const int *const numTargs,
1678  const QMPCVec grpPntChrgVec,
1679  QMPCVec &grpAppldChrgVec) {
1680 
1681  DebugM(1,"Processing QM-MM bonds in rank zero.\n");
1682 
1683  // Indices and IDs of MM atoms which will be passed as point charges.
1684  std::map<int, int> mmPntChrgMap ;
1685 
1686  // Build the map of atom ID and charge for this QM group's point charges.
1687  // They may be changed by different charge distribution schemes and/or
1688  // dummy aotm placements.
1689  for (size_t i=0; i<grpPntChrgVec.size(); i++) {
1690 
1691  mmPntChrgMap.insert(std::pair<int,int>(grpPntChrgVec[i].id, (int) i) );
1692 
1693  grpAppldChrgVec.push_back( grpPntChrgVec[i] ) ;
1694 
1695  }
1696 
1697  // If we are treating QM-MM bonds and if there are any
1698  // in this QM group, we have to treat the MM partial charge
1699  // using some scheme.
1700  for (int bondIt = 0; bondIt < numBonds; bondIt++) {
1701 
1702  // Gets the global index of this particular QM-MM bond.
1703  int bondIndx = qmGrpBondsGrp[bondIt] ;
1704 
1705  auto retIt = mmPntChrgMap.find(qmMMBondedIndxGrp[bondIt]) ;
1706 
1707  // Checks if this MM atom was included as a
1708  // point charge due proximity to a QM atom.
1709  if (retIt == mmPntChrgMap.end()) {
1710  // If it wasn't, there is an error somwhere.
1711 
1712  iout << iERROR << "The MM atom " << qmMMBondedIndxGrp[bondIt]
1713  << " is bound to a QM atom, but it was not selected as a poitn charge."
1714  << " Check your cutoff radius!\n" << endi ;
1715 
1716  NAMD_die("Charge placement error in QM-MM bond.");
1717  }
1718 
1719  // Gets the (local) index of the MM atom in the point charge vector.
1720  int mmIndex = (*retIt).second;
1721  // Gets the position of the MM atom and its partial charge.
1722  Position mmPos = grpAppldChrgVec[mmIndex].position ;
1723  BigReal mmCharge = grpAppldChrgVec[mmIndex].charge/numTargs[bondIndx] ;
1724 
1725  // gives part of the MM charge to neighboring atoms
1726  for (int i=0; i<numTargs[bondIndx]; i++){
1727 
1728  int targetIndxGLobal = chargeTarget[bondIndx][i] ;
1729 
1730  retIt = mmPntChrgMap.find(targetIndxGLobal);
1731 
1732  // If one of the neighboring atoms which should receive charge
1733  // is not among partial charges, stop and run away.
1734  if (retIt == mmPntChrgMap.end()) {
1735 
1736  iout << iERROR << "The MM atom " << targetIndxGLobal
1737  << " is bound to the MM atom of a QM-MM bond and is needed for"
1738  << " the required bond scheme"
1739  << " but it was not selected as a poitn charge."
1740  << " Check your cutoff radius!\n" << endi ;
1741 
1742  NAMD_die("Charge placement error in QM-MM bond.");
1743  }
1744 
1745  int trgIndxLocal = (*retIt).second;
1746 
1747  // Choose charge treatment method and apply it.
1748  switch (simParams->qmBondScheme) {
1749 
1750  // Charge Schiftig scheme.
1751  case QMSCHEMECS:
1752  {
1753 
1754  // Charge Shifting Scheme (DOI: 10.1021/ct100530r)
1755  // Here we create a dipole to counter act the charge movement
1756  // we created by moving parts of the MM charge to target MM atoms.
1757  //
1758  // The MM1 charge is equally divided and added to all MM2 atoms.
1759  // Two point charges are created in the direction of the MM1-MM2 bond,
1760  // one before and one after the MM2 atom.
1761  //
1762  // Below we see the diagram of the positions of new charges along
1763  // the direction of the bond between the MM1 atom and a target
1764  // MM2 atom, wiht respect to the new Dummy atom (a Hydrogen).
1765  //
1766  // QM -------- MM1(p0) ------------ MM2
1767  // QM ------ H ------------ (+p0) - (MM2 +p0) - (-p0)
1768 
1769  Position trgPos = grpAppldChrgVec[trgIndxLocal].position ;
1770 
1771  // We now store in MM2 to which PC it is bound (the index for MM1).
1772  // This will be used to decompose virtual PC forces onto MM1 and MM2.
1773  grpAppldChrgVec[trgIndxLocal].qmGrpID = -2;
1774  grpAppldChrgVec[trgIndxLocal].homeIndx = mmIndex;
1775 
1776  // We create a new point charge at the same position so that
1777  // the fraction of charge from the MM1 atom is "added" to the
1778  // MM2 position, but without actually changing the charge of
1779  // the original MM2 atom. This trick helps keeping the original
1780  // charge for electrostatic calculations after the new charges
1781  // of QM atoms is received from the QM software.
1782  Real Cq0 = 1.0;
1783  grpAppldChrgVec.push_back(
1784  // Cq is stored in the distance slot, since it has not been computed
1785  // for a virtual charge.
1786  ComputeQMPntChrg(trgPos, mmCharge, trgIndxLocal, 0, -1, Cq0, 0, 0)
1787  );
1788 
1789  Vector bondVec = trgPos - mmPos ;
1790 
1791  // For Cq-plus virtual charge
1792  Real Cqp = 0.94;
1793  // For Cq-minus virtual charge
1794  Real Cqm = 1.06;
1795  Vector bondVec1 = bondVec*Cqp ;
1796  Vector bondVec2 = bondVec*Cqm ;
1797 
1798  Position chrgPos1 = mmPos + bondVec1;
1799  Position chrgPos2 = mmPos + bondVec2;
1800 
1801  BigReal trgChrg1 = mmCharge;
1802  BigReal trgChrg2 = -1*mmCharge;
1803 
1804  grpAppldChrgVec.push_back(
1805  ComputeQMPntChrg(chrgPos1, trgChrg1, trgIndxLocal, 0, -1, Cqp, 0, 0)
1806  );
1807 
1808  grpAppldChrgVec.push_back(
1809  ComputeQMPntChrg(chrgPos2, trgChrg2, trgIndxLocal, 0, -1, Cqm, 0, 0)
1810  );
1811 
1812  } break;
1813 
1814  // Redistributed Charge and Dipole scheme.
1815  case QMSCHEMERCD:
1816  {
1817 
1818 // grpAppldChrgVec[trgIndxLocal].charge -= mmCharge ;
1819 
1820  // Redistributed Charge and Dipole method (DOI: 10.1021/ct100530r)
1821  // Here we create a dipole to counter act the charge movement
1822  // we created by moving parts of the MM charge to target MM atoms.
1823  //
1824  // The MM1 charge is equally divided and subtracted from all MM2 atoms.
1825  // One point charge is created in the midpoint of the MM1-MM2 bond.
1826  //
1827  // Below we see the diagram of the positions of new charges along
1828  // the direction of the bond between the MM1 atom and a target
1829  // MM2 atom, wiht respect to the new Dummy atom (a Hydrogen).
1830  //
1831  // QM -------- MM1(p0) -------------- MM2
1832  // QM ------ H ------------ (+2*p0) - (MM2 -p0)
1833 
1834  Position trgPos = grpAppldChrgVec[trgIndxLocal].position ;
1835 
1836  // We now store in MM2 to which PC it is bound (the index for MM1).
1837  // This will be used to decompose virtual PC forces onto MM1 and MM2.
1838  grpAppldChrgVec[trgIndxLocal].qmGrpID = -2;
1839  grpAppldChrgVec[trgIndxLocal].homeIndx = mmIndex;
1840 
1841  // We create a new point charge at the same position so that
1842  // the fraction of charge from the MM1 atom is "added" to the
1843  // MM2 position, but without actually changing the charge of
1844  // the original MM2 atom. This trick helps keeping the original
1845  // charge for electrostatic calculations after the new charges
1846  // of QM atoms is received from the QM software.
1847  Real Cq0 = 1.0;
1848  grpAppldChrgVec.push_back(
1849  // Cq is stored in the distance slot, since it has not been computed
1850  // for a virtual charge.
1851  ComputeQMPntChrg(trgPos, -1*mmCharge, trgIndxLocal, 0, -1, Cq0, 0, 0)
1852  );
1853 
1854  Vector bondVec = trgPos - mmPos ;
1855 
1856  // For Cq-plus virtual charge
1857  Real Cq1 = 0.5;
1858  Vector bondVec1 = bondVec*Cq1 ;
1859 
1860  Position chrgPos1 = mmPos + bondVec1;
1861 
1862  BigReal trgChrg1 = 2*mmCharge;
1863 
1864  grpAppldChrgVec.push_back(
1865  ComputeQMPntChrg(chrgPos1, trgChrg1, trgIndxLocal, 0, -1, Cq1, 0, 0)
1866  );
1867 
1868 
1869 
1870  } break;
1871 
1872  // The following schemes simply make the charges of the
1873  // MM atom(s) equal to zero. No redistribution or
1874  // artificial dipole is attmepted.
1875  //
1876  // For Z1, only the MM1 atom has its charge set to zero.
1877  // For Z2, the MM1 and all MM2 atom charges are changed.
1878  // For Z3, the MM1, all MM2 and all MM3 atom charges are changed.
1879  //
1880  // All modification are local ONLY. They do not change the atom's
1881  // charge for the classical side of the simulation. Only the QM
1882  // side ceases to see the electrostatic influence of these atoms,
1883  // and they cease to see the QM region electrostatics as well.
1884 
1885  // Z1, Z2 and Z3 schemes do the same thing. The only difference
1886  // is the target list, which is created in MoleculeQM.C.
1887  case QMSCHEMEZ1:
1888 
1889  case QMSCHEMEZ2:
1890 
1891  case QMSCHEMEZ3:
1892  grpAppldChrgVec[trgIndxLocal].charge = 0.0;
1893  break;
1894  }
1895 
1896  }
1897 
1898  // We keep this "point charge" so we can calculate forces on it later
1899  // but it will not influence the QM system.
1900  // We use the qmGrpID variable to send the message that this point charge
1901  // should be ignored since this variable will not be relevant anymore.
1902  // All point charges gathered here are for a specific qmGroup.
1903  grpAppldChrgVec[mmIndex].qmGrpID = -1 ;
1904  }
1905 
1906  return ;
1907 }
1908 
1909 
1911 
1912  // All the preparation that used to be in this function was moved to
1913  // recvPartQM, which is called first in rank zero.
1914 
1915  if (noPC) {
1916  // Even if we have QM-MM bonds, the point charge messages
1917  // are handled in recvPartQM
1918  delete msg;
1919  }
1920  else {
1921  pntChrgCoordMsgs[numArrivedPntChrgMsg] = msg;
1922  ++numArrivedPntChrgMsg;
1923 
1924  if ( numArrivedPntChrgMsg < numSources )
1925  return;
1926  }
1927 
1928  // Resets counters for next step.
1929  numRecQMRes = 0;
1930 
1931  totalEnergy = 0;
1932 
1933  for ( int k=0; k<3; ++k )
1934  for ( int l=0; l<3; ++l )
1935  totVirial[k][l] = 0;
1936 
1937  // ALL DATA ARRIVED --- CALCULATE FORCES
1938 
1939  const char *const *const elementArray = molPtr->get_qmElements() ;
1940  const char *const *const dummyElmArray = molPtr->get_qmDummyElement();
1941  const int *const qmGrpSize = molPtr->get_qmGrpSizes();
1942 
1943  const BigReal *const dummyBondVal = molPtr->get_qmDummyBondVal();
1944  const int *const grpNumBonds = molPtr->get_qmGrpNumBonds() ;
1945  const int *const *const qmMMBond = molPtr->get_qmMMBond() ;
1946  const int *const *const qmGrpBonds = molPtr->get_qmGrpBonds() ;
1947  const int *const *const qmMMBondedIndx = molPtr->get_qmMMBondedIndx() ;
1948  const int *const *const chargeTarget = molPtr->get_qmMMChargeTarget() ;
1949  const int *const numTargs = molPtr->get_qmMMNumTargs() ;
1950 
1951  BigReal constants = COULOMB * simParams->nonbondedScaling / (simParams->dielectric) ;
1952  // COULOMB is in kcal*Angs/(mol*e^2)
1953 // BigReal constants = COULOMB ;
1954 
1955  if ( qmPCFreq > 0 ) {
1956  DebugM(4,"Using point charge stride of " << qmPCFreq << "\n")
1957  // In case we are skiping steps between PC selection, store a main list
1958  // in rank zero for future steps. Then we only update positions untill
1959  // we reach a new "qmPCFreq" step, when a new PC selection is made.
1960 
1961  if (timeStep % qmPCFreq == 0) {
1962  DebugM(4,"Loading a new selection of PCs.\n")
1963 
1964  // We only re-set this arrya in a step where a new PC selection is made.
1965  pntChrgMsgVec.clear();
1966  for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
1967  // Accumulates the message point charges into a local vector.
1968  for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
1969  pntChrgMsgVec.push_back(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
1970  }
1971  }
1972 
1973  // This fast array is created to send the entire PC IDs list to all
1974  // PEs in the next step.
1975  pcIDListSize = pntChrgMsgVec.size();
1976  pcIDList = new int[pcIDListSize] ;
1977  for (size_t i=0; i<pcIDListSize; i++) {
1978  pcIDList[i] = pntChrgMsgVec[i].id;
1979 
1980  // Loads home indexes of MM atoms received as point charges.
1981  // Each force receives the home index of its atom with respect to the
1982  // local set of atoms in each node.
1983  force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
1984  }
1985  }
1986  else {
1987  DebugM(4,"Updating position/homeIdex of old PC selection.\n")
1988 
1989  // We build a sorted array so that we can quickly access the unordered
1990  // data we just received, and update positions of the same selection
1991  // of point charges.
1992  pcDataSort.clear();
1993  for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
1994  // Accumulates the message point charges into a local sorted array.
1995  for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
1996  pcDataSort.load(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
1997  }
1998  }
1999  pcDataSort.sort();
2000 
2001  iout << "Loaded new positions in sorted array: " << pcDataSort.size() << "\n" << endi;
2002 
2003  // If we are using a set of point charges that was selected in a
2004  // previous step, we update the PC POSITION ONLY.
2005  for (size_t i=0; i<pntChrgMsgVec.size(); i++) {
2006 
2007  auto pntr = pcDataSort.find(pntChrgMsgVec[i]);
2008 
2009  pntChrgMsgVec[i].position = pntr->position ;
2010  pntChrgMsgVec[i].homeIndx = pntr->homeIndx ;
2011 
2012  // Loads home indexes of MM atoms received as point charges.
2013  // Each force receives the home index of its atom with respect to the
2014  // local set of atoms in each node.
2015  force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
2016  }
2017  }
2018  }
2019  else {
2020  DebugM(4,"Updating PCs at every step.\n")
2021 
2022  pntChrgMsgVec.clear();
2023  for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
2024  // Accumulates the message point charges into a local vector.
2025  for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
2026  pntChrgMsgVec.push_back(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
2027  }
2028  }
2029 
2030  // Loads home indexes of MM atoms received as point charges.
2031  for (size_t i=0; i<pntChrgMsgVec.size(); i++) {
2032  // Each force receives the home index of its atom with respect to the
2033  // local set of atoms in each node.
2034 
2035  #ifdef DEBUG_QM
2036  if (force[pntChrgMsgVec[i].id].homeIndx != -1
2037  and force[pntChrgMsgVec[i].id].homeIndx != pntChrgMsgVec[i].homeIndx
2038  ) {
2039  iout << iERROR << "Overloading point charge "
2040  << pntChrgMsgVec[i].id << "; home index: "
2041  << force[pntChrgMsgVec[i].id].homeIndx << " with " << pntChrgMsgVec[i].homeIndx
2042  << "\n" << endi ;
2043  NAMD_die("Error preparing QM simulations.");
2044  }
2045  #endif
2046 
2047  force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
2048  }
2049  }
2050 
2051  // Reset counter for next timestep
2052  numArrivedPntChrgMsg = 0;
2053 
2054  DebugM(4,"A total of " << pntChrgMsgVec.size()
2055  << " point charges arrived." << std::endl);
2056 
2057  DebugM(4,"Starting QM groups processing.\n");
2058 
2059  QMAtmVec grpQMAtmVec;
2060  QMPCVec grpPntChrgVec;
2061 
2062  // Final set of charges, created or not, that is sent to the QM software.
2063  // This set will depend on how QM-MM bonds are processed and presented to the
2064  // QM region.
2065  QMPCVec grpAppldChrgVec;
2066 
2067  // Vector of dummy atoms created to treat QM-MM bonds.
2068  std::vector<dummyData> dummyAtoms ;
2069 
2070  // Initializes the loop for receiving the QM results.
2071  thisProxy[0].recvQMResLoop() ;
2072 
2073  // Iterator for target PE where QM simulations will run.
2074  int peIter = 0;
2075 
2076  for (int grpIter = 0; grpIter < numQMGrps; grpIter++) {
2077 
2078  grpQMAtmVec.clear();
2079  grpPntChrgVec.clear();
2080  grpAppldChrgVec.clear();
2081  dummyAtoms.clear();
2082 
2083  DebugM(4,"Calculating QM group " << grpIter +1
2084  << " (ID: " << grpID[grpIter] << ")." << std::endl);
2085 
2086  DebugM(4,"Compiling Config Lines into one string for message...\n");
2087 
2088  // This will hold a big sting with all configuration lines the user supplied.
2089  int lineIter = 0 ;
2090  std::string configLines ;
2091  StringList *current = Node::Object()->configList->find("QMConfigLine");
2092  for ( ; current; current = current->next ) {
2093  std::string confLineStr(current->data);
2094 
2095  // In case we need to add charges to MOPAC command line.
2096  if (simParams->qmFormat == QMFormatMOPAC && simParams->qmMOPACAddConfigChrg && lineIter == 0) {
2097  std::ostringstream itosConv ;
2098  itosConv << grpChrg[grpIter] ;
2099  confLineStr.append( " CHARGE=" );
2100  confLineStr.append( itosConv.str() );
2101 
2102  }
2103 
2104  configLines.append(confLineStr);
2105  configLines.append("\n");
2106 
2107  lineIter++;
2108  }
2109 
2110  DebugM(4,"Determining point charges...\n");
2111 
2112  Real qmTotalCharge = 0;
2113  // Loads the QM atoms for this group.
2114  for (int qmIt=0; qmIt<numQMAtms; qmIt++){
2115  if (qmCoord[qmIt].qmGrpID == grpID[grpIter]) {
2116  grpQMAtmVec.push_back(qmCoord[qmIt]);
2117  qmTotalCharge += qmCoord[qmIt].charge;
2118  }
2119  }
2120  if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
2121  qmTotalCharge = roundf(qmTotalCharge) ;
2122  }
2123 
2124  // Sorts the vector so that the QM message is always built with the same order of atoms.
2125  std::sort(grpQMAtmVec.begin(), grpQMAtmVec.end(), custom_ComputeQMAtom_Less);
2126 
2127  Real pcTotalCharge = 0;
2128  // Loads the point charges to a local vector for this QM group.
2129  for (auto pntChrgIt = pntChrgMsgVec.begin();
2130  pntChrgIt != pntChrgMsgVec.end(); pntChrgIt++) {
2131  if ((*pntChrgIt).qmGrpID == grpID[grpIter] ) {
2132  grpPntChrgVec.push_back( (*pntChrgIt) );
2133  pcTotalCharge += (*pntChrgIt).charge;
2134  }
2135  }
2136  if ((fabsf(roundf(pcTotalCharge) - pcTotalCharge) <= 0.001f) ) {
2137  pcTotalCharge = roundf(pcTotalCharge) ;
2138  }
2139 
2140  #ifdef DEBUG_QM
2141  if (grpQMAtmVec.size() != qmGrpSize[grpIter]) {
2142  iout << iERROR << "The number of QM atoms received for group " << grpID[grpIter]
2143  << " does not match the expected: " << grpQMAtmVec.size()
2144  << " vs " << qmGrpSize[grpIter] << "\n" << endi ;
2145 
2146  NAMD_die("Error processing QM group.");
2147  }
2148  #endif
2149 
2150  DebugM(1,"Found " << grpPntChrgVec.size() << " point charges for QM group "
2151  << grpIter << " (ID: " << grpID[grpIter] << "; Num QM atoms: "
2152  << grpQMAtmVec.size() << "; Num QM-MM bonds: "
2153  << grpNumBonds[grpIter] << ")" << std::endl);
2154 
2155  DebugM(1,"Total QM charge: " << qmTotalCharge
2156  << "; Total point-charge charge: " << pcTotalCharge << std::endl);
2157 
2158  // If we have a frequency for LSS update, check if we shoudl do it in
2159  // the current time step.
2160  if ( qmLSSFreq > 0 && ((timeStep + 1) % qmLSSFreq == 0 )) {
2161  lssUpdate(grpIter, grpQMAtmVec, grpPntChrgVec);
2162  }
2163 
2164  // This function checks data and treats the charge (and existence) of
2165  // the MM atoms in and around QM-MM bonds. It is only executed in
2166  // electrostatic embeding QM/MM simulations.
2167  if (! noPC ) {
2168  procBonds(grpNumBonds[grpIter], qmGrpBonds[grpIter],
2169  qmMMBondedIndx[grpIter],
2170  chargeTarget, numTargs,
2171  grpPntChrgVec, grpAppldChrgVec) ;
2172  }
2173  else {
2174  grpAppldChrgVec = grpPntChrgVec;
2175  }
2176 
2177  // For future use, we get the pairs of indexes of QM atoms and MM atoms which are
2178  // bound in QM-MM bonds.
2179  std::vector< std::pair<int,int> > qmPCLocalIndxPairs ;
2180 
2181  // Create and position dummy atoms.
2182  Position mmPos(0), qmPos(0);
2183  for (int dummyIt = 0; dummyIt < grpNumBonds[grpIter]; dummyIt++) {
2184 
2185  int qmMMBondIndx = qmGrpBonds[grpIter][dummyIt] ;
2186 
2187  BigReal bondVal = dummyBondVal[qmMMBondIndx] ;
2188 
2189  int mmAtmIndx = qmMMBond[qmMMBondIndx][0] ;
2190  int qmAtmIndx = qmMMBond[qmMMBondIndx][1] ;
2191 
2192  // Sicne we don't know in which patch/node the QM atom is, or the
2193  // order in which they will arrive in rank zero, we have
2194  // no direct index to it.
2195  #ifdef DEBUG_QM
2196  bool missingQM = true, missingMM = true;
2197  #endif
2198  size_t qmIt ;
2199  for (qmIt=0; qmIt<grpQMAtmVec.size(); qmIt++){
2200  if (grpQMAtmVec[qmIt].id == qmAtmIndx) {
2201  qmPos = grpQMAtmVec[qmIt].position;
2202 
2203  #ifdef DEBUG_QM
2204  missingQM = false;
2205  #endif
2206 
2207  break;
2208  }
2209  }
2210  // The same is true about the MM atom to which the QM atom is bound,
2211  // we must look
2212  size_t pcIt;
2213  for (pcIt=0; pcIt < grpPntChrgVec.size(); pcIt++) {
2214  if (grpPntChrgVec[pcIt].id == mmAtmIndx) {
2215  mmPos = grpPntChrgVec[pcIt].position ;
2216 
2217  #ifdef DEBUG_QM
2218  missingMM = false;
2219  #endif
2220 
2221  break;
2222  }
2223  }
2224 
2225  qmPCLocalIndxPairs.push_back(std::pair<int,int>(qmIt, pcIt) );
2226 
2227  #ifdef DEBUG_QM
2228  // Checks if the MM atom was included as a
2229  // point charge due proximity to a QM atom, and if the QM atom arrived.
2230  if ( missingMM or missingQM ) {
2231  // If it wasn't, there is an error somwhere.
2232 
2233  if (missingMM) {
2234  iout << iERROR << "The MM atom " << mmAtmIndx
2235  << " is bound to a QM atom, but it was not selected as a poitn charge."
2236  << " Check your cutoff radius!\n" << endi ;
2237 
2238  NAMD_die("Error in QM-MM bond processing.");
2239  }
2240  if (missingQM) {
2241  iout << iERROR << "The QM atom " << qmAtmIndx
2242  << " is bound to an MM atom, but it was not sent to rank zero for processing."
2243  << " Check your configuration!\n" << endi ;
2244 
2245  NAMD_die("Error in QM-MM bond processing.");
2246  }
2247  }
2248  #endif
2249 
2250  Vector bondVec = mmPos - qmPos ;
2251 
2252  if (bondValType == QMLENTYPE) {
2253  // If a length is defined by the user, or a default len
2254  // is used, we calculate the unit vector for the displacement
2255  // and multiply by the desired length in order
2256  // to get the final dummy atom position relative to the
2257  // QM atom.
2258  bondVec = bondVec.unit() ;
2259  bondVec *= bondVal ;
2260  }
2261  else if (bondValType == QMRATIOTYPE) {
2262  // If distance a ratio was defined by the user, then
2263  // the displacement vector is multiplied by that ratio
2264  // to get the final dummy atom position relative to the
2265  // QM atom.
2266  bondVec *= bondVal ;
2267  }
2268 
2269  Position dummyPos = qmPos + bondVec;
2270 
2271  DebugM(1,"Creating dummy atom " << dummyPos << " ; QM ind: "
2272  << qmIt << " ; PC ind: " << pcIt << std::endl);
2273 
2274  dummyAtoms.push_back(dummyData(dummyPos, qmIt, qmMMBondIndx)) ;
2275 
2276  }
2277 
2278  DebugM(3, "Creating data for " << grpQMAtmVec.size() << " QM atoms "
2279  << dummyAtoms.size() << " dummy atoms " << grpPntChrgVec.size()
2280  << " real point charges and " << grpAppldChrgVec.size() - grpPntChrgVec.size()
2281  << " virtual point charges\n") ;
2282 
2283  int dataSize = grpQMAtmVec.size() + dummyAtoms.size() + grpAppldChrgVec.size();
2284  QMGrpCalcMsg *msg = new (dataSize, configLines.size(), 0) QMGrpCalcMsg;
2285  msg->timestep = timeStep;
2286  msg->grpIndx = grpIter;
2287  msg->peIter = peIter;
2288  msg->charge = grpChrg[grpIter];
2289  msg->multiplicity = grpMult[grpIter];
2290  msg->numQMAtoms = grpQMAtmVec.size();
2291  msg->numAllAtoms = grpQMAtmVec.size() + dummyAtoms.size();
2292  msg->numRealPntChrgs = grpPntChrgVec.size(); // The original set of classical atoms.
2293  msg->numAllPntChrgs = grpAppldChrgVec.size(); // The extended set with virtual point charges.
2294  msg->secProcOn = simParams->qmSecProcOn ;
2295  msg->constants = constants;
2296  msg->PMEOn = simParams->PMEOn ;
2297  if (msg->PMEOn)
2298  msg->PMEEwaldCoefficient = simParams->PMEEwaldCoefficient ;
2299  msg->switching = simParams->qmPCSwitchOn;
2300  msg->switchType = simParams->qmPCSwitchType;
2301  msg->cutoff = simParams->cutoff;
2302  msg->swdist = simParams->switchingDist;
2303  msg->pcScheme = simParams->qmPCScheme;
2304  msg->qmAtmChrgMode = simParams->qmChrgMode;
2305 
2306  strncpy(msg->baseDir, simParams->qmBaseDir, 256);
2307  strncpy(msg->execPath, simParams->qmExecPath, 256);
2308  if (msg->secProcOn)
2309  strncpy(msg->secProc, simParams->qmSecProc, 256);
2310 
2311  if (simParams->qmPrepProcOn && (timeStep == simParams->firstTimestep)) {
2312  msg->prepProcOn = true;
2313  strncpy(msg->prepProc, simParams->qmPrepProc, 256);
2314  } else
2315  msg->prepProcOn = false;
2316 
2317  QMAtomData *dataP = msg->data;
2318 
2319  for (int i=0; i<grpQMAtmVec.size(); i++) {
2320  dataP->position = grpQMAtmVec[i].position ;
2321  dataP->charge = grpQMAtmVec[i].charge ;
2322  dataP->id = grpQMAtmVec[i].id ;
2323  dataP->bountToIndx = -1;
2324  dataP->type = QMATOMTYPE_QM ;
2325  strncpy(dataP->element,elementArray[grpQMAtmVec[i].id],3);
2326  dataP++;
2327  }
2328 
2329  for (int i=0; i<dummyAtoms.size(); i++) {
2330  dataP->position = dummyAtoms[i].pos ;
2331  dataP->charge = 0 ;
2332  dataP->id = -1 ;
2333  dataP->bountToIndx = dummyAtoms[i].qmInd ;
2334  dataP->type = QMATOMTYPE_DUMMY ;
2335  strncpy(dataP->element,dummyElmArray[dummyAtoms[i].bondIndx],3);
2336  dataP++;
2337  }
2338 
2339  for (int i=0; i<grpAppldChrgVec.size(); i++) {
2340  dataP->position = grpAppldChrgVec[i].position ;
2341  dataP->charge = grpAppldChrgVec[i].charge ;
2342  // Point charges created to handle QM-MM bonds will have an id of -1.
2343  dataP->id = grpAppldChrgVec[i].id ;
2344  dataP->bountToIndx = -1;
2345  dataP->dist = grpAppldChrgVec[i].dist ;
2346  // If we are loading the classical atoms' charges
2347  // the point charge type is 1, unless it is from an
2348  // atom which is bound to a QM atom.
2349  if (i < grpPntChrgVec.size()) {
2350  if (grpAppldChrgVec[i].qmGrpID == -1) {
2351  dataP->type = QMPCTYPE_IGNORE ;
2352  }
2353  else if (grpAppldChrgVec[i].qmGrpID == -2) {
2354  dataP->type = QMPCTYPE_CLASSICAL ;
2355  dataP->bountToIndx = grpAppldChrgVec[i].homeIndx;
2356  }
2357  else {
2358  dataP->type = QMPCTYPE_CLASSICAL ;
2359  }
2360  }
2361  else {
2362  // Extra charges are created to handle QM-MM bonds (if they exist).
2363  dataP->type = QMPCTYPE_EXTRA ;
2364  dataP->bountToIndx = grpAppldChrgVec[i].id;
2365  }
2366  dataP++;
2367  }
2368 
2369  QMAtomData *qmP = msg->data ;
2370  QMAtomData *pcP = msg->data + msg->numAllAtoms ;
2371 
2372  // With this, every QM atom knows to which MM atom is is bound,
2373  // and vice versa. This will be usefull later on to prevent them from
2374  // feeling eachother's electrostatic charges AND to place the dummy
2375  // atom forces on the "real" atoms that form the bond.
2376  for( auto vecPtr = qmPCLocalIndxPairs.begin();
2377  vecPtr != qmPCLocalIndxPairs.end();
2378  vecPtr++ ) {
2379 
2380  int qmLocInd = (*vecPtr).first;
2381  int pcLocInd = (*vecPtr).second;
2382 
2383  qmP[qmLocInd].bountToIndx = pcLocInd ;
2384  pcP[pcLocInd].bountToIndx = qmLocInd ;
2385  }
2386 
2387 
2388  strcpy(msg->configLines, configLines.c_str());
2389 
2390  if (cSMDon)
2391  calcCSMD(grpIter, msg->numQMAtoms, qmP, cSMDForces) ;
2392 
2393  int targetPE = qmPEs[peIter] ;
2394 
2395  DebugM(4,"Sending QM group " << grpIter << " (ID " << grpID[grpIter]
2396  << ") to PE " << targetPE << std::endl);
2397 
2398  switch (simParams->qmFormat) {
2399  // Creates the command line that will be executed according to the
2400  // chosen QM software, as well as the input file with coordinates.
2401  case QMFormatORCA:
2402  QMProxy[targetPE].calcORCA(msg) ;
2403  break;
2404 
2405  case QMFormatMOPAC:
2406  QMProxy[targetPE].calcMOPAC(msg) ;
2407  break;
2408 
2409  case QMFormatUSR:
2410  QMProxy[targetPE].calcUSR(msg) ;
2411  break;
2412  }
2413 
2414  peIter++;
2415 
2416  if (peIter == qmPEs.size())
2417  peIter = 0;
2418  }
2419 
2420 }
2421 
2423 
2424 // iout << iINFO << "Storing QM results for region " << resMsg->grpIndx
2425 // << " (ID: " << grpID[resMsg->grpIndx]
2426 // << ") with original energy: " << endi;
2427 // std::cout << std::fixed << std::setprecision(6) << resMsg->energyOrig << endi;
2428 // iout << " Kcal/mol\n" << endi;
2429 
2430 // if (resMsg->energyCorr != resMsg->energyOrig) {
2431 // iout << iINFO << "PME corrected energy: " << endi;
2432 // std::cout << std::fixed << std::setprecision(6) << resMsg->energyCorr << endi;
2433 // iout << " Kcal/mol\n" << endi;
2434 // }
2435 
2436  if ( (timeStep % simParams->qmEnergyOutFreq ) == 0 ) {
2437  iout << QMETITLE(timeStep);
2438  iout << FORMAT(grpID[resMsg->grpIndx]);
2439  iout << FORMAT(resMsg->energyOrig);
2440  if (resMsg->energyCorr != resMsg->energyOrig) iout << FORMAT(resMsg->energyCorr);
2441  iout << "\n\n" << endi;
2442  }
2443 
2444  totalEnergy += resMsg->energyCorr ;
2445 
2446  for ( int k=0; k<3; ++k )
2447  for ( int l=0; l<3; ++l )
2448  totVirial[k][l] += resMsg->virial[k][l];
2449 
2450  QMForce *fres = resMsg->force ;
2451  Real qmTotalCharge = 0;
2452 
2453  for (int i=0; i<resMsg->numForces; i++) {
2454 
2455  force[ fres[i].id ].force += fres[i].force;
2456 
2457  // Indicates the result is a QM atom, and we should get its updated charge.
2458  if (fres[i].replace == 1) {
2459  force[ fres[i].id ].charge = fres[i].charge;
2460  qmTotalCharge += fres[i].charge;
2461  }
2462  }
2463 
2464  if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
2465  qmTotalCharge = roundf(qmTotalCharge) ;
2466  }
2467 
2468  // In case we are calculating cSMD forces, apply them now.
2469  if (cSMDon) {
2470  for ( int i=0; i < cSMDindxLen[resMsg->grpIndx]; i++ ) {
2471  int cSMDit = cSMDindex[resMsg->grpIndx][i];
2472  force[ cSMDpairs[cSMDit][0] ].force += cSMDForces[cSMDit];
2473  }
2474  }
2475 
2476  DebugM(4,"QM total charge received is " << qmTotalCharge << std::endl);
2477 
2478  DebugM(4,"Current accumulated energy is " << totalEnergy << std::endl);
2479 
2480  numRecQMRes++;
2481 
2482  delete resMsg;
2483 }
2484 
2486 
2487  // Writes a DCD file with the charges of all QM atoms at a frequency
2488  // defined by the user in qmOutFreq.
2489  if ( simParams->qmOutFreq > 0 &&
2490  timeStep % simParams->qmOutFreq == 0 ) {
2491 
2492  iout << iINFO << "Writing QM charge output at step "
2493  << timeStep << "\n" << endi;
2494 
2495  Real *x = outputData,
2496  *y = outputData + molPtr->get_numQMAtoms(),
2497  *z = outputData + 2*molPtr->get_numQMAtoms();
2498 
2499  for (int qmIt=0; qmIt<numQMAtms; qmIt++){
2500  x[qmIt] = qmCoord[qmIt].id;
2501  y[qmIt] = force[qmCoord[qmIt].id].charge ;
2502  z[qmIt] = 0;
2503  }
2504 
2505  write_dcdstep(dcdOutFile, numQMAtms, x, y, z, 0) ;
2506  }
2507 
2508  // Writes a DCD file with the positions of all QM atoms at a frequency
2509  // defined by the user in qmPosOutFreq.
2510  if ( simParams->qmPosOutFreq > 0 &&
2511  timeStep % simParams->qmPosOutFreq == 0 ) {
2512 
2513  iout << iINFO << "Writing QM position output at step "
2514  << timeStep << "\n" << endi;
2515 
2516  SortedArray<idIndxStr> idIndx;
2517 
2518  for(int i=0; i<numQMAtms;i++) {
2519  idIndx.insert( idIndxStr(qmCoord[i].id, i) );
2520  }
2521  idIndx.sort();
2522 
2523  Real *x = outputData,
2524  *y = outputData + molPtr->get_numQMAtoms(),
2525  *z = outputData + 2*molPtr->get_numQMAtoms();
2526 
2527  for (int qmIt=0; qmIt<numQMAtms; qmIt++){
2528  x[qmIt] = qmCoord[idIndx[qmIt].indx].position.x;
2529  y[qmIt] = qmCoord[idIndx[qmIt].indx].position.y;
2530  z[qmIt] = qmCoord[idIndx[qmIt].indx].position.z;
2531  }
2532 
2533  write_dcdstep(dcdPosOutFile, numQMAtms, x, y, z, 0) ;
2534  }
2535 
2536  // distribute forces
2537  DebugM(4,"Distributing QM forces for all ranks.\n");
2538  for ( int j=0; j < numSources; ++j ) {
2539 
2540  std::set<int> auxset0;
2541  DebugM(1,"Sending forces and charges to source " << j << std::endl);
2542 
2543  QMCoordMsg *qmmsg = 0;
2544  QMPntChrgMsg *pcmsg = 0 ;
2545 
2546  int totForces = 0;
2547  int sourceNode = -1;
2548 
2549  if (pntChrgCoordMsgs == NULL) {
2550 
2551  qmmsg = qmCoordMsgs[j];
2552  qmCoordMsgs[j] = 0;
2553 
2554  totForces = qmmsg->numAtoms ;
2555 
2556  sourceNode = qmmsg->sourceNode;
2557  }
2558  else {
2559  pcmsg = pntChrgCoordMsgs[j];
2560  pntChrgCoordMsgs[j] = 0;
2561 
2562  sourceNode = pcmsg->sourceNode;
2563 
2564  // Since we receive two different messages from nodes, there is no
2565  // guarantee the two sets of messages will come in the same order.
2566  // Therefore, we match the messages by comaring their sourceNodes.
2567  for (int aux=0; aux<numSources; aux++) {
2568 
2569  if (qmCoordMsgs[aux] == 0)
2570  continue;
2571 
2572  qmmsg = qmCoordMsgs[aux];
2573 
2574  if (qmmsg->sourceNode == sourceNode) {
2575  qmCoordMsgs[aux] = 0;
2576  break;
2577  }
2578  }
2579 
2580  DebugM(1,"Building force mesage for rank "
2581  << pcmsg->sourceNode << std::endl);
2582 
2583  totForces = qmmsg->numAtoms + pcmsg->numAtoms;
2584 
2585  // I want to count number of different atom ids
2586  // which is the size of force array.
2587  // Avoids double counting of forces
2588  for ( int i=0; i < qmmsg->numAtoms; ++i ) {
2589  auxset0.insert(qmmsg->coord[i].id);
2590  }
2591  for ( int i=0; i < pcmsg->numAtoms; ++i ) {
2592  auxset0.insert(pcmsg->coord[i].id);
2593  }
2594  totForces = auxset0.size(); // Fixes same patch different QM regions double counting
2595  }
2596 
2597  QMForceMsg *fmsg = new (totForces, subsArray.size(), 0) QMForceMsg;
2598  fmsg->PMEOn = simParams->PMEOn;
2599  fmsg->numForces = totForces;
2600  fmsg->numLssDat = subsArray.size();
2601 
2602  DebugM(1,"Loading QM forces.\n");
2603 
2604  // This iterator is used in BOTH loops!
2605  int forceIter = 0;
2606  auxset0.clear(); // Need to keep track of forces by id that have been copied to fmsg
2607  for ( int i=0; i < qmmsg->numAtoms; ++i ) {
2608  fmsg->force[forceIter] = force[qmmsg->coord[i].id];
2609  forceIter++;
2610  auxset0.insert(qmmsg->coord[i].id); // This should add each qm atom once
2611  }
2612 
2613  delete qmmsg;
2614 
2615  if (pntChrgCoordMsgs != NULL) {
2616  DebugM(1,"Loading PntChrg forces.\n");
2617 
2618  for ( int i=0; i < pcmsg->numAtoms; ++i ) {
2619  // (QM atoms that are PC or repeated PC) are not copied to fmsg
2620  if (auxset0.find(pcmsg->coord[i].id) == auxset0.end()) {
2621  fmsg->force[forceIter] = force[pcmsg->coord[i].id];
2622  forceIter++;
2623  auxset0.insert(pcmsg->coord[i].id);
2624  }
2625  }
2626 
2627  delete pcmsg;
2628  }
2629 
2630  DebugM(1,"A total of " << forceIter << " forces were loaded." << std::endl);
2631 
2632  for ( int i=0; i < subsArray.size(); ++i ) {
2633  fmsg->lssDat[i] = subsArray[i];
2634  }
2635 
2636  #ifdef DEBUG_QM
2637  if (subsArray.size() > 0)
2638  DebugM(3,"A total of " << subsArray.size() << " LSS data structures were loaded." << std::endl);
2639  #endif
2640 
2641  if ( ! j ) {
2642  fmsg->energy = totalEnergy;
2643  for ( int k=0; k<3; ++k )
2644  for ( int l=0; l<3; ++l )
2645  fmsg->virial[k][l] = totVirial[k][l];
2646  } else {
2647  fmsg->energy = 0;
2648  for ( int k=0; k<3; ++k )
2649  for ( int l=0; l<3; ++l )
2650  fmsg->virial[k][l] = 0;
2651  }
2652 
2653  DebugM(4,"Sending forces...\n");
2654 
2655  QMProxy[sourceNode].recvForce(fmsg);
2656 
2657  }
2658 
2659  DebugM(4,"All forces sent from node zero.\n");
2660 }
2661 
2663 
2664  if (CkMyPe()) {
2665  for (int i=0; i<fmsg->numLssDat; i++) {
2666  subsArray.add(fmsg->lssDat[i]) ;
2667  }
2668  }
2669 
2670  QMCompute->saveResults(fmsg);
2671 }
2672 
2674 
2676 
2677  bool callReplaceForces = false;
2678 
2679  int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
2680  int numQMGrps = Node::Object()->molecule->get_qmNumGrps();
2681  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
2682  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
2683  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
2684 
2685  int totalNumbAtoms = 0;
2686  for (ap = ap.begin(); ap != ap.end(); ap++) {
2687  totalNumbAtoms += (*ap).p->getNumAtoms();
2688  }
2689 
2690  // This is kept to be deleted in the next step so that the message can be
2691  // used in "replaceForces" routine later on, in case there are forces to
2692  // be replaced. The "replaceForces" function uses but DOES NOT free the pointer,
2693  // so we free the data from the previous iteration and allocate a new one for
2694  // the current iteration (remember that the number of atoms can change in a
2695  // home patch between iterations).
2696  if (oldForces != 0)
2697  delete [] oldForces;
2698  oldForces = new ExtForce[totalNumbAtoms] ;
2699 
2700  for (int i=0; i < totalNumbAtoms; ++i) {
2701  oldForces[i].force = Force(0) ;
2702  }
2703 
2704  DebugM(1,"Force array has been created and zeroed in rank "
2705  << CkMyPe() << std::endl);
2706 
2707  DebugM(1,"Preparing " << fmsg->numForces << " forces in rank "
2708  << CkMyPe() << std::endl);
2709 
2710  QMForce *results_ptr = fmsg->force;
2711  // Iterates over the received forces and place them on the full array.
2712  for (int i=0; i < fmsg->numForces; ++i, ++results_ptr) {
2713  // For now we may have more than one item in the message acting on the same
2714  // atom, such as an MM atom which is a point charge for two or more QM regions.
2715 
2716  oldForces[results_ptr->homeIndx].force += results_ptr->force;
2717  oldForces[results_ptr->homeIndx].replace = results_ptr->replace;
2718 
2719  if (results_ptr->replace == 1)
2720  callReplaceForces = true;
2721 
2722  // If the atom is in a QM group, update its charge to the local (this homePatch)
2723  // copy of the qmAtmChrg array.
2724  if (qmAtomGroup[results_ptr->id] > 0 && (fmsg->PMEOn || (numQMGrps > 1) ) ) {
2725 
2726  // Loops over all QM atoms (in all QM groups) comparing their global indices
2727  for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
2728 
2729  if (qmAtmIndx[qmIter] == results_ptr->id) {
2730  qmAtmChrg[qmIter] = results_ptr->charge;
2731  break;
2732  }
2733 
2734  }
2735 
2736  }
2737 
2738  }
2739 
2740  DebugM(1,"Placing forces on force boxes in rank "
2741  << CkMyPe() << std::endl);
2742 
2743  // Places the received forces on the force array for each patch.
2744  int homeIndxIter = 0;
2745  for (ap = ap.begin(); ap != ap.end(); ap++) {
2746  Results *r = (*ap).forceBox->open();
2747  Force *f = r->f[Results::normal];
2748  const CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
2749  int localNumAtoms = (*ap).p->getNumAtoms();
2750 
2751  for(int i=0; i<localNumAtoms; ++i) {
2752 
2753  f[i] += oldForces[homeIndxIter].force;
2754 
2755  ++homeIndxIter;
2756  }
2757 
2758  if ( callReplaceForces )
2759  (*ap).p->replaceForces(oldForces);
2760 
2761  (*ap).forceBox->close(&r);
2762 
2763  }
2764 
2765  DebugM(1,"Forces placed on force boxes in rank "
2766  << CkMyPe() << std::endl);
2767 
2768  if (fmsg->PMEOn) {
2769 
2770  DebugM(1,"PME ON! Accessing PMEmgr in rank " << CkMyPe() << std::endl);
2771 
2772  ComputePmeMgr *mgrP = CProxy_ComputePmeMgr::ckLocalBranch(
2773  CkpvAccess(BOCclass_group).computePmeMgr) ;
2774 
2775  DebugM(4, "Initiating ComputePme::doQMWork on rank " << CkMyPe() << " over "
2776  << getComputes(mgrP).size() << " pmeComputes." << std::endl) ;
2777 
2778  for ( int i=0; i< getComputes(mgrP).size(); ++i ) {
2779  // WorkDistrib::messageEnqueueWork(pmeComputes[i]);
2780  getComputes(mgrP)[i]->doQMWork();
2781  }
2782  }
2783 
2784  DebugM(1,"Submitting reduction in rank " << CkMyPe() << std::endl);
2785 
2786  reduction->item(REDUCTION_MISC_ENERGY) += fmsg->energy;
2787  reduction->item(REDUCTION_VIRIAL_NORMAL_XX) += fmsg->virial[0][0];
2788  reduction->item(REDUCTION_VIRIAL_NORMAL_XY) += fmsg->virial[0][1];
2789  reduction->item(REDUCTION_VIRIAL_NORMAL_XZ) += fmsg->virial[0][2];
2790  reduction->item(REDUCTION_VIRIAL_NORMAL_YX) += fmsg->virial[1][0];
2791  reduction->item(REDUCTION_VIRIAL_NORMAL_YY) += fmsg->virial[1][1];
2792  reduction->item(REDUCTION_VIRIAL_NORMAL_YZ) += fmsg->virial[1][2];
2793  reduction->item(REDUCTION_VIRIAL_NORMAL_ZX) += fmsg->virial[2][0];
2794  reduction->item(REDUCTION_VIRIAL_NORMAL_ZY) += fmsg->virial[2][1];
2795  reduction->item(REDUCTION_VIRIAL_NORMAL_ZZ) += fmsg->virial[2][2];
2796  reduction->submit();
2797 
2798  delete fmsg ;
2799 }
2800 
2802 {
2803  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
2804  FILE *inputFile,*outputFile,*chrgFile;
2805  int iret;
2806 
2807  const size_t lineLen = 256;
2808  char *line = new char[lineLen];
2809 
2810  std::string qmCommand, inputFileName, outputFileName, pntChrgFileName;
2811 
2812  // For coulomb calculation
2813  BigReal constants = msg->constants;
2814 
2815  double gradient[3];
2816 
2817  DebugM(4,"Running MOPAC on PE " << CkMyPe() << std::endl);
2818 
2819  QMAtomData *atmP = msg->data ;
2820  QMAtomData *pcP = msg->data + msg->numAllAtoms ;
2821 
2822  // We create a pointer for PME correction, which may point to
2823  // a copy of the original point charge array, with unchanged charges,
2824  // or points to the original array in case no switching or charge
2825  // scheme is being used.
2826  QMAtomData *pcPpme = NULL;
2827  if (msg->switching) {
2828 
2829  if (msg->PMEOn)
2830  pcPpme = new QMAtomData[msg->numRealPntChrgs];
2831 
2832  pntChrgSwitching(msg, pcPpme) ;
2833  } else {
2834  pcPpme = pcP;
2835  }
2836 
2837  int retVal = 0;
2838  struct stat info;
2839 
2840  // For each QM group, create a subdirectory where files will be palced.
2841  std::string baseDir(msg->baseDir);
2842  std::ostringstream itosConv ;
2843  if ( CmiNumPartitions() > 1 ) {
2844  baseDir.append("/") ;
2845  itosConv << CmiMyPartition() ;
2846  baseDir += itosConv.str() ;
2847  itosConv.str("");
2848  itosConv.clear() ;
2849 
2850  if (stat(msg->baseDir, &info) != 0 ) {
2851  CkPrintf( "Node %d cannot access directory %s\n",
2852  CkMyPe(), baseDir.c_str() );
2853  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
2854  }
2855  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
2856  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
2857  retVal = mkdir(baseDir.c_str(), S_IRWXU);
2858  }
2859 
2860  }
2861  baseDir.append("/") ;
2862  itosConv << msg->grpIndx ;
2863  baseDir += itosConv.str() ;
2864 
2865  if (stat(msg->baseDir, &info) != 0 ) {
2866  CkPrintf( "Node %d cannot access directory %s\n",
2867  CkMyPe(), baseDir.c_str() );
2868  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
2869  }
2870  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
2871  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
2872  retVal = mkdir(baseDir.c_str(), S_IRWXU);
2873  }
2874 
2875  #ifdef DEBUG_QM
2876  Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
2877  #endif
2878 
2879  inputFileName.clear();
2880  inputFileName.append(baseDir.c_str()) ;
2881  inputFileName.append("/qmmm_") ;
2882  inputFileName += itosConv.str() ;
2883  inputFileName.append(".input") ;
2884 
2885  // Opens file for coordinate and parameter input
2886  inputFile = fopen(inputFileName.c_str(),"w");
2887  if ( ! inputFile ) {
2888  iout << iERROR << "Could not open input file for writing: "
2889  << inputFileName << "\n" << endi ;
2890  NAMD_err("Error writing QM input file.");
2891  }
2892 
2893  // Builds the command that will be executed
2894  qmCommand.clear();
2895  qmCommand.append("cd ");
2896  qmCommand.append(baseDir);
2897  qmCommand.append(" ; ");
2898  qmCommand.append(msg->execPath) ;
2899  qmCommand.append(" ") ;
2900  qmCommand.append(inputFileName) ;
2901  qmCommand.append(" > /dev/null 2>&1") ;
2902 
2903  // Builds the file name where MOPAC will place the gradient
2904  // This is also relative to the input file name.
2905  outputFileName = inputFileName ;
2906  outputFileName.append(".aux") ;
2907 
2908  if (msg->numAllPntChrgs) {
2909  // Builds the file name where we will write the point charges.
2910  pntChrgFileName.clear();
2911  pntChrgFileName.append(baseDir.c_str()) ;
2912  pntChrgFileName.append("/mol.in") ;
2913 
2914  chrgFile = fopen(pntChrgFileName.c_str(),"w");
2915  if ( ! chrgFile ) {
2916  iout << iERROR << "Could not open charge file for writing: "
2917  << pntChrgFileName << "\n" << endi ;
2918  NAMD_err("Error writing point charge file.");
2919  }
2920 
2921  iret = fprintf(chrgFile,"\n%d %d\n",
2922  msg->numQMAtoms, msg->numAllAtoms - msg->numQMAtoms);
2923  if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
2924  }
2925 
2926  // writes configuration lines to the input file.
2927  std::stringstream ss(msg->configLines) ;
2928  std::string confLineStr;
2929  while (std::getline(ss, confLineStr) ) {
2930  confLineStr.append("\n");
2931  iret = fprintf(inputFile,confLineStr.c_str());
2932  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
2933  }
2934 
2935 
2936  iret = fprintf(inputFile,"\n");
2937  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
2938 
2939  DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file "
2940  << inputFileName.c_str() << " and " << msg->numAllPntChrgs
2941  << " point charges in file " << pntChrgFileName.c_str() << "\n");
2942 
2943  // write QM and dummy atom coordinates to input file and
2944  // MM electric field from MM point charges.
2945  for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
2946 
2947  double x = atmP->position.x;
2948  double y = atmP->position.y;
2949  double z = atmP->position.z;
2950 
2951  iret = fprintf(inputFile,"%s %f 1 %f 1 %f 1\n",
2952  atmP->element,x,y,z);
2953  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
2954 
2955  if (msg->numAllPntChrgs) {
2956  BigReal phi = 0;
2957 
2958  // The Electrostatic Potential is calculated according to
2959  // the QMMM Keyword documentation for MOPAC
2960  // http://openmopac.net/manual/QMMM.html
2961  pcP = msg->data + msg->numAllAtoms ;
2962  for ( size_t j=0; j < msg->numAllPntChrgs; ++j, ++pcP ) {
2963 
2964  // In case of QM-MM bonds, the charge of the MM1 atom is ignored
2965  if (pcP->type == QMPCTYPE_IGNORE)
2966  continue;
2967 
2968  double charge = pcP->charge;
2969 
2970  double xMM = pcP->position.x;
2971  double yMM = pcP->position.y;
2972  double zMM = pcP->position.z;
2973 
2974  BigReal rij = sqrt((x-xMM)*(x-xMM) +
2975  (y-yMM)*(y-yMM) +
2976  (z-zMM)*(z-zMM) ) ;
2977 
2978  phi += charge/rij ;
2979  }
2980 
2981  // We use the same Coulomb constant used in the rest of NAMD
2982  // instead of the suggested rounded 332 suggested by MOPAC.
2983  phi = phi*constants ;
2984 
2985  iret = fprintf(chrgFile,"%s %f %f %f %f\n",
2986  atmP->element,x,y,z, phi);
2987  if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
2988  }
2989  }
2990 
2991  DebugM(4,"Closing input and charge file\n");
2992 
2993  if (msg->numAllPntChrgs)
2994  fclose(chrgFile);
2995 
2996  fclose(inputFile);
2997 
2998  if (msg->prepProcOn) {
2999 
3000  std::string prepProc(msg->prepProc) ;
3001  prepProc.append(" ") ;
3002  prepProc.append(inputFileName) ;
3003  iret = system(prepProc.c_str());
3004  if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
3005  if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
3006  }
3007 
3008  // runs QM command
3009  DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
3010  iret = system(qmCommand.c_str());
3011 
3012  if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
3013  if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
3014 
3015  if (msg->secProcOn) {
3016 
3017  std::string secProc(msg->secProc) ;
3018  secProc.append(" ") ;
3019  secProc.append(inputFileName) ;
3020  itosConv.str("");
3021  itosConv.clear() ;
3022  itosConv << msg->timestep ;
3023  secProc.append(" ") ;
3024  secProc += itosConv.str() ;
3025 
3026  iret = system(secProc.c_str());
3027  if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
3028  if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
3029  }
3030 
3031  // remove coordinate file
3032 // iret = remove(inputFileName);
3033 // if ( iret ) { NAMD_err(0); }
3034 
3035  // remove coordinate file
3036 // iret = remove(pntChrgFileName);
3037 // if ( iret ) { NAMD_err(0); }
3038 
3039  // opens output file
3040  DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
3041  outputFile = fopen(outputFileName.c_str(),"r");
3042  if ( ! outputFile ) {
3043  iout << iERROR << "Could not find QM output file!\n" << endi;
3044  NAMD_err(0);
3045  }
3046 
3047  // Resets the pointers.
3048  atmP = msg->data ;
3049  pcP = msg->data + msg->numAllAtoms ;
3050 
3051  // Allocates the return message.
3052  QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
3053  resMsg->grpIndx = msg->grpIndx;
3054  resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
3055  resMsg->energyOrig = 0;
3056  resMsg->energyCorr = 0;
3057  for ( int k=0; k<3; ++k )
3058  for ( int l=0; l<3; ++l )
3059  resMsg->virial[k][l] = 0;
3060  QMForce *resForce = resMsg->force ;
3061 
3062  // We split the data into two pointers so we can "skip" the dummy atoms
3063  // which have no representation in NAMD.
3064  for (int i=0; i<resMsg->numForces; i++) {
3065  resForce[i].force = 0;
3066  resForce[i].charge = 0 ;
3067  if (i < msg->numQMAtoms) {
3068  // We use the replace field to indicate QM atoms,
3069  // which will have charge information.
3070  resForce[i].replace = 1 ;
3071  resForce[i].id = atmP->id;
3072  atmP++;
3073  }
3074  else {
3075  // We use the replace field to indicate QM atoms,
3076  // which will have charge information.
3077  resForce[i].replace = 0 ;
3078  resForce[i].id = pcP->id;
3079  pcP++;
3080  }
3081  }
3082 
3083  // Resets the pointers.
3084  atmP = msg->data ;
3085  pcP = msg->data + msg->numAllAtoms ;
3086 
3087  // Reads the data form the output file created by the QM software.
3088  // Gradients over the QM atoms, and Charges for QM atoms will be read.
3089 
3090  size_t atmIndx = 0, gradCount = 0;
3091  Bool gradFields = false, chargeFields = false;
3092  Bool chargesRead = false, gradsRead = false;
3093  while ( fgets(line, lineLen, outputFile) != NULL){
3094  // We loop over the lines of the output file untill we find
3095  // the line that initiates the "atom charges" lines. Then
3096  // we read all lines and expect to get one or more charges
3097  // per line, separated by space(s), untill we find a line that
3098  // initiates the "gradients", then once more, we expect several
3099  // numbers separated by space(s). When the "overlap matrix"
3100  // string is found, we break the loop and stop reading the file.
3101 
3102  // if ( strstr(line,"TOTAL_ENERGY") != NULL ) {
3103  if ( strstr(line,"HEAT_OF_FORMATION") != NULL ) {
3104 
3105  char strEnergy[14], *endPtr ;
3106 
3107  //strncpy(strEnergy, line + 17, 13) ;
3108  strncpy(strEnergy, line + 28, 13) ;
3109  strEnergy[13] = '\0';
3110 
3111  // We have to convert the notation from FORTRAN double precision to
3112  // the natural, normal, reasonable and not terribly out dated format.
3113  resMsg->energyOrig = strtod(strEnergy, &endPtr);
3114  if ( *endPtr == 'D' ) {
3115  *endPtr = 'e';
3116  resMsg->energyOrig = strtod(strEnergy, &endPtr);
3117  }
3118 
3119  // In MOPAC, the total energy is given in EV, so we convert to Kcal/mol
3120 // resMsg->energyOrig *= 23.060 ; // We now read Heat of Formation, which is given in Kcal/Mol
3121 
3122 // DebugM(4,"Reading QM energy from file: " << resMsg->energyOrig << "\n");
3123 
3124  resMsg->energyCorr = resMsg->energyOrig;
3125 
3126  continue;
3127  }
3128 
3129  if ( strstr(line,"ATOM_CHARGES") != NULL ) {
3130  gradFields = false;
3131  chargeFields = true;
3132  atmIndx = 0;
3133 
3134  // If the user asked for charges NOT to be read from MOPAC,
3135  // we skip the charge reading step.
3136  if (msg->qmAtmChrgMode == QMCHRGNONE) {
3137  chargeFields = false;
3138  atmIndx = msg->numAllAtoms;
3139  }
3140  else {
3141  chargeFields = true;
3142  atmIndx = 0;
3143  }
3144 
3145  // Now we expect the following line(s) to have atom charges
3146  continue;
3147  }
3148 
3149  if ( strstr(line,"GRADIENTS") != NULL ) {
3150 
3151  // Now that all charges have been read, checks if the
3152  // numbers match
3153  if (atmIndx != msg->numAllAtoms) {
3154  NAMD_die("Error reading QM forces file. Wrong number of atom charges");
3155  }
3156 
3157  chargesRead = true;
3158 
3159  // Now we expect the following line(s) to have gradients
3160  chargeFields = false ;
3161  gradFields = true;
3162  gradCount = 0;
3163  atmIndx = 0;
3164 
3165  continue;
3166  }
3167 
3168  if ( strstr(line,"OVERLAP_MATRIX") != NULL ) {
3169 
3170  // Now that all gradients have been read, checks if the
3171  // numbers match
3172  if (atmIndx != msg->numAllAtoms) {
3173  NAMD_die("Error reading QM forces file. Wrong number of gradients");
3174  }
3175 
3176  gradsRead = true;
3177 
3178  // Nothing more to read from the ".aux" file
3179  break;
3180  }
3181 
3182  char result[20] ;
3183  size_t strIndx = 0;
3184 
3185  if (chargeFields) {
3186  while ((strIndx < (strlen(line)-9)) && (strlen(line)-1 >=9 ) ) {
3187 
3188  strncpy(result, line+strIndx,9) ;
3189  result[9] = '\0';
3190 
3191  Real localCharge = atof(result);
3192 
3193  // If we are reading charges from QM atoms, store them.
3194  if (atmIndx < msg->numQMAtoms ) {
3195  atmP[atmIndx].charge = localCharge;
3196  resForce[atmIndx].charge = localCharge;
3197  }
3198 
3199  // If we are reading charges from Dummy atoms,
3200  // place them on the appropriate QM atoms.
3201  if ( msg->numQMAtoms <= atmIndx &&
3202  atmIndx < msg->numAllAtoms ) {
3203  // We keep the calculated charge in the dummy atom
3204  // so that we can calculate electrostatic forces later on.
3205  atmP[atmIndx].charge = localCharge;
3206 
3207  // The dummy atom points to the QM atom to which it is bound.
3208  int qmInd = atmP[atmIndx].bountToIndx ;
3209  resForce[qmInd].charge += localCharge;
3210  }
3211 
3212  strIndx += 9;
3213  atmIndx++;
3214 
3215  // If we found all charges for QM and dummy atoms, break the loop
3216  // and stop reading the line.
3217  if (atmIndx == msg->numAllAtoms) {
3218  chargeFields = false;
3219  break;
3220  }
3221 
3222  }
3223 
3224  }
3225 
3226  int gradLength ; // Change for variable length MOPAC output
3227  if (gradFields) {
3228  if (atmIndx == 0) {
3229  double buf[10];
3230  int numfirstline = sscanf(line,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
3231  &buf[0],&buf[1],&buf[2],&buf[3],&buf[4],&buf[5],
3232  &buf[6],&buf[7],&buf[8],&buf[9]);
3233  gradLength = strlen(line)/numfirstline;
3234  }
3235  while ((strIndx < (strlen(line)-gradLength)) && (strlen(line)-1 >= gradLength ) ) {
3236 
3237  strncpy(result, line+strIndx,gradLength) ;
3238  result[gradLength] = '\0';
3239 
3240  gradient[gradCount] = atof(result);
3241  if (gradCount == 2) {
3242 
3243  if (atmIndx < msg->numQMAtoms ) {
3244  // Gradients in MOPAC are written in kcal/mol/A.
3245  resForce[atmIndx].force.x = -1*gradient[0];
3246  resForce[atmIndx].force.y = -1*gradient[1];
3247  resForce[atmIndx].force.z = -1*gradient[2];
3248  }
3249 
3250  // If we are reading forces applied on Dummy atoms,
3251  // place them on the appropriate QM and MM atoms to conserve energy.
3252 
3253  // This implementation was based on the description in
3254  // DOI: 10.1002/jcc.20857 : Equations 30 to 32
3255  if ( msg->numQMAtoms <= atmIndx &&
3256  atmIndx < msg->numAllAtoms ) {
3257  // The dummy atom points to the QM atom to which it is bound.
3258  // The QM atom and the MM atom (in a QM-MM bond) point to each other.
3259  int qmInd = atmP[atmIndx].bountToIndx ;
3260  int mmInd = atmP[qmInd].bountToIndx ;
3261 
3262  Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
3263  Real mmqmDist = dir.length() ;
3264 
3265  Real linkDist = Vector(atmP[atmIndx].position -
3266  atmP[qmInd].position).length() ;
3267 
3268  Force mmForce(0), qmForce(0),
3269  linkForce(gradient[0], gradient[1], gradient[2]);
3270  linkForce *= -1;
3271 
3272  Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
3273  // Unit vectors
3274  Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
3275  Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
3276  Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
3277  Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
3278 
3279  qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv +
3280  (xDelta)*base) )*xuv;
3281 
3282  qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv +
3283  (yDelta)*base) )*yuv;
3284 
3285  qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv +
3286  (zDelta)*base) )*zuv;
3287 
3288 
3289  mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
3290  (xDelta)*base) )*xuv;
3291 
3292  mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
3293  (yDelta)*base) )*yuv;
3294 
3295  mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
3296  (zDelta)*base) )*zuv;
3297 
3298  resForce[qmInd].force += qmForce;
3299  resForce[msg->numQMAtoms + mmInd].force += mmForce;
3300  }
3301 
3302  gradCount = 0;
3303  atmIndx++;
3304  } else {
3305  gradCount++;
3306  }
3307 
3308  strIndx += gradLength;
3309 
3310  // If we found all gradients for QM atoms, break the loop
3311  // and keep the next gradient line from being read, if there
3312  // is one. Following gradients, if present, will refer to link
3313  // atoms, and who cares about those?.
3314  if (atmIndx == msg->numAllAtoms) {
3315  gradFields = false;
3316  break;
3317  }
3318  }
3319  }
3320 
3321  }
3322 
3323  delete [] line;
3324 
3325  fclose(outputFile);
3326 
3327  // In case charges are not to be read form the QM software output,
3328  // we load the origianl atom charges.
3329  if (msg->qmAtmChrgMode == QMCHRGNONE) {
3330  int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
3331  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
3332  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
3333 
3334  atmIndx = 0 ;
3335  for (; atmIndx < msg->numQMAtoms; atmIndx++) {
3336 
3337  // Loops over all QM atoms (in all QM groups) comparing their global indices
3338  for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
3339 
3340  if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
3341 
3342  atmP[atmIndx].charge = qmAtmChrg[qmIter];
3343  resForce[atmIndx].charge = qmAtmChrg[qmIter];
3344 
3345  break;
3346  }
3347 
3348  }
3349 
3350  }
3351  for (size_t j=msg->numQMAtoms; j<msg->numAllAtoms; ++j ) {
3352  atmP[j].charge = 0;
3353  }
3354  }
3355 
3356  // remove force file
3357 // DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
3358 // iret = remove(outputFileName);
3359 // if ( iret ) { NAMD_err(0); }
3360 
3361  if (! (chargesRead && gradsRead) ) {
3362  NAMD_die("Error reading QM forces file. Not all data could be read!");
3363  }
3364 
3365  DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
3366 
3367  atmP = msg->data ;
3368  pcP = msg->data + msg->numAllAtoms ;
3369 
3370  // The initial point charge index for the *force* message is the number of QM
3371  // atoms, since the dummy atoms have no representation in NAMD
3372  int pcIndx = msg->numQMAtoms;
3373 
3374  // ---- WARNING ----
3375  // We need to apply the force felt by each QM atom due to the classical
3376  // point charges sent to MOPAC.
3377  // MOPAC takes the MM electrostatic potential into cosideration ONLY
3378  // when performing the SCF calculation and when returning the energy
3379  // of the system, but it does not apply the potential to the final
3380  // gradient calculation, therefore, we calculate the resulting force
3381  // here.
3382  // Initialize force vector for electrostatic contribution
3383  std::vector<Force> qmElecForce ;
3384  for (size_t j=0; j<msg->numAllAtoms; ++j )
3385  qmElecForce.push_back( Force(0) );
3386 
3387  // We loop over point charges and distribute the forces applied over
3388  // virtual point charges to the MM1 and MM2 atoms (if any virtual PCs exists).
3389  for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
3390 
3391  // No force was applied to the QM region due to this charge, so it
3392  // does not receive any back from the QM region. It must be an MM1 atom
3393  // from a QM-MM bond.
3394  if (pcP[i].type == QMPCTYPE_IGNORE)
3395  continue;
3396 
3397  Force totalForce(0);
3398 
3399  BigReal pntCharge = pcP[i].charge;
3400 
3401  Position posMM = pcP[i].position ;
3402 
3403  for (size_t j=0; j<msg->numAllAtoms; ++j ) {
3404 
3405  BigReal qmCharge = atmP[j].charge ;
3406 
3407  BigReal force = pntCharge*qmCharge*constants ;
3408 
3409  Position rVec = posMM - atmP[j].position ;
3410 
3411  force /= rVec.length2();
3412 
3413  // We accumulate the total force felt by a point charge
3414  // due to the QM charge distribution. This total force
3415  // will be applied to the point charge if it is a real one,
3416  // or will be distirbuted over MM1 and MM2 point charges, it
3417  // this is a virtual point charge.
3418  totalForce += force*rVec.unit();
3419 
3420  // Accumulate force felt by a QM atom due to point charges.
3421  qmElecForce[j] += -1*force*rVec.unit();
3422  }
3423 
3424  if (pcP[i].type == QMPCTYPE_CLASSICAL) {
3425  // Checking pcP was not a QM atom in another region
3426  // If so the interaction is accounted in the other region
3427  if (qmAtomGroup[pcP[i].id] == 0) {
3428  // If we already ignored MM1 charges, we take all other
3429  // non-virtual charges and apply forces directly to them.
3430  resForce[pcIndx].force += totalForce;
3431  }
3432  }
3433  else {
3434  // If we are handling virtual PC, we distribute the force over
3435  // MM1 and MM2.
3436 
3437  // Virtual PC are bound to MM2.
3438  int mm2Indx = pcP[i].bountToIndx;
3439  // MM2 charges are bound to MM1.
3440  int mm1Indx = pcP[mm2Indx].bountToIndx;
3441 
3442  Real Cq = pcP[i].dist;
3443 
3444  Force mm1Force = (1-Cq)*totalForce ;
3445  Force mm2Force = Cq*totalForce ;
3446 
3447  resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
3448  resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
3449  }
3450 
3451  }
3452 
3453  // We now apply the accumulated electrostatic force contribution
3454  // to QM atoms.
3455  for (size_t j=0; j<msg->numAllAtoms; ++j ) {
3456 
3457  if (j < msg->numQMAtoms ) {
3458  resForce[j].force += qmElecForce[j];
3459 
3460  } else {
3461  // In case the QM atom is a Link atom, we redistribute
3462  // its force as bove.
3463 
3464  // The dummy atom points to the QM atom to which it is bound.
3465  // The QM atom and the MM atom (in a QM-MM bond) point to each other.
3466  int qmInd = atmP[j].bountToIndx ;
3467  int mmInd = atmP[qmInd].bountToIndx ;
3468 
3469  Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
3470  Real mmqmDist = dir.length() ;
3471 
3472  Real linkDist = Vector(atmP[j].position -
3473  atmP[qmInd].position).length() ;
3474 
3475  Force linkForce = qmElecForce[j];
3476 
3477  Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
3478  // Unit vectors
3479  Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
3480  Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
3481  Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
3482  Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
3483 
3484  Force qmForce = (linkForce*((1 - linkDist/mmqmDist)*xuv +
3485  (xDelta)*base) )*xuv;
3486 
3487  qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv +
3488  (yDelta)*base) )*yuv;
3489 
3490  qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv +
3491  (zDelta)*base) )*zuv;
3492 
3493 
3494  Force mmForce = (linkForce*((linkDist/mmqmDist)*xuv -
3495  (xDelta)*base) )*xuv;
3496 
3497  mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
3498  (yDelta)*base) )*yuv;
3499 
3500  mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
3501  (zDelta)*base) )*zuv;
3502 
3503  resForce[qmInd].force += qmForce;
3504  resForce[msg->numQMAtoms + mmInd].force += mmForce;
3505  }
3506  }
3507 
3508  // Adjusts forces from PME, canceling contributions from the QM and
3509  // direct Coulomb forces calculated here.
3510  if (msg->PMEOn) {
3511 
3512  DebugM(1,"Correcting forces and energy for PME.\n");
3513 
3514  ewaldcof = msg->PMEEwaldCoefficient;
3515  BigReal TwoBySqrtPi = 1.12837916709551;
3516  pi_ewaldcof = TwoBySqrtPi * ewaldcof;
3517 
3518  for (size_t i=0; i < msg->numQMAtoms; i++) {
3519 
3520  BigReal p_i_charge = atmP[i].charge ;
3521  Position pos_i = atmP[i].position ;
3522 
3523  const BigReal kq_i = p_i_charge * constants;
3524 
3525  for (size_t j=i+1; j < msg->numQMAtoms; j++) {
3526 
3527  BigReal p_j_charge = atmP[j].charge ;
3528 
3529  Position pos_j = atmP[j].position ;
3530 
3531  BigReal r = Vector(pos_i - pos_j).length() ;
3532 
3533  BigReal tmp_a = r * ewaldcof;
3534  BigReal tmp_b = erfc(tmp_a);
3535  BigReal corr_energy = tmp_b;
3536  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
3537 
3538 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
3539  BigReal recip_energy = (1-tmp_b)/r;
3540 
3541 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
3542  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
3543 
3544  // Final force and energy correction for this pair of atoms.
3545  BigReal energy = kq_i * p_j_charge * recip_energy ;
3546 
3547  Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
3548 
3549  // The force is *subtracted* from the total force acting on
3550  // both atoms. The sign on fixForce corrects the orientation
3551  // of the subtracted force.
3552 // DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
3553 // << std::endl);
3554 // DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
3555 // << std::endl);
3556 // DebugM(4,"Force correction: " << fixForce << std::endl);
3557 // DebugM(4,"Energy correction: " << energy << "\n");
3558  resForce[i].force -= fixForce ;
3559  resForce[j].force -= -1*fixForce ;
3560 
3561  // The energy is *subtracted* from the total energy calculated here.
3562  resMsg->energyCorr -= energy;
3563  }
3564  }
3565 
3566 // DebugM(4,"Corrected energy QM-QM interactions: " << resMsg->energyCorr << "\n");
3567 
3568  pcIndx = msg->numQMAtoms;
3569  // We only loop over point charges from real atoms, ignoring the ones
3570  // created to handle QM-MM bonds.
3571  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
3572 
3573  BigReal p_i_charge = pcPpme[i].charge;
3574  Position pos_i = pcPpme[i].position ;
3575 
3576  const BigReal kq_i = p_i_charge * constants;
3577 
3578  Force fixForce = 0;
3579 
3580  for (size_t j=0; j<msg->numQMAtoms; ++j ) {
3581 
3582  BigReal p_j_charge = atmP[j].charge ;
3583 
3584  Position pos_j = atmP[j].position ;
3585 
3586  BigReal r = Vector(pos_i - pos_j).length() ;
3587 
3588  BigReal tmp_a = r * ewaldcof;
3589  BigReal tmp_b = erfc(tmp_a);
3590  BigReal corr_energy = tmp_b;
3591  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
3592 
3593 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
3594  BigReal recip_energy = (1-tmp_b)/r;
3595 
3596 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
3597  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
3598 
3599  // Final force and energy correction for this pair of atoms.
3600  BigReal energy = kq_i * p_j_charge * recip_energy ;
3601 
3602  fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
3603 
3604  // The energy is *subtracted* from the total energy calculated here.
3605  resMsg->energyCorr -= energy;
3606  }
3607 
3608  // The force is *subtracted* from the total force acting on
3609  // the point charge..
3610 // DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
3611 // << std::endl);
3612 // DebugM(4,"Force correction: " << fixForce << std::endl);
3613  resForce[pcIndx].force -= kq_i*fixForce ;
3614  }
3615 
3616  }
3617 
3618  // Calculates the virial contribution form the forces on QM atoms and
3619  // point charges calculated here.
3620  for (size_t i=0; i < msg->numQMAtoms; i++) {
3621  // virial used by NAMD is -'ve of normal convention, so reverse it!
3622  // virial[i][j] in file should be sum of -1 * f_i * r_j
3623  for ( int k=0; k<3; ++k )
3624  for ( int l=0; l<3; ++l )
3625  resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
3626  }
3627 
3628  pcIndx = msg->numQMAtoms; // Index in the real PC force array.
3629  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
3630  // virial used by NAMD is -'ve of normal convention, so reverse it!
3631  // virial[i][j] in file should be sum of -1 * f_i * r_j
3632  for ( int k=0; k<3; ++k )
3633  for ( int l=0; l<3; ++l )
3634  resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
3635  }
3636 
3637 
3638 
3639  // Send message to rank zero with results.
3640  QMProxy[0].recvQMRes(resMsg);
3641 
3642  if (msg->switching && msg->PMEOn)
3643  delete [] pcPpme;
3644 
3645  delete msg;
3646  return ;
3647 }
3648 
3650 {
3651  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
3652  FILE *inputFile,*outputFile,*chrgFile;
3653  int iret;
3654 
3655  const size_t lineLen = 256;
3656  char *line = new char[lineLen];
3657 
3658  std::string qmCommand, inputFileName, outputFileName, pntChrgFileName, pcGradFilename;
3659  std::string tmpRedirectFileName, pcGradFileName;
3660 
3661  // For coulomb calculation
3662  BigReal constants = msg->constants;
3663 
3664  double gradient[3];
3665 
3666  DebugM(4,"Running ORCA on PE " << CkMyPe() << std::endl);
3667 
3668  QMAtomData *pcP = msg->data + msg->numAllAtoms ;
3669 
3670  // We create a pointer for PME correction, which may point to
3671  // a copy of the original point charge array, with unchanged charges,
3672  // or points to the original array in case no switching or charge
3673  // scheme is being used.
3674  QMAtomData *pcPpme = NULL;
3675  if (msg->switching) {
3676 
3677  if (msg->PMEOn)
3678  pcPpme = new QMAtomData[msg->numRealPntChrgs];
3679 
3680  pntChrgSwitching(msg, pcPpme) ;
3681  } else {
3682  pcPpme = pcP;
3683  }
3684 
3685  int retVal = 0;
3686  struct stat info;
3687 
3688  // For each QM group, create a subdirectory where files will be palced.
3689  std::string baseDir(msg->baseDir);
3690  std::ostringstream itosConv ;
3691  if ( CmiNumPartitions() > 1 ) {
3692  baseDir.append("/") ;
3693  itosConv << CmiMyPartition() ;
3694  baseDir += itosConv.str() ;
3695  itosConv.str("");
3696  itosConv.clear() ;
3697 
3698  if (stat(msg->baseDir, &info) != 0 ) {
3699  CkPrintf( "Node %d cannot access directory %s\n",
3700  CkMyPe(), baseDir.c_str() );
3701  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
3702  }
3703  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
3704  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
3705  retVal = mkdir(baseDir.c_str(), S_IRWXU);
3706  }
3707 
3708  }
3709  baseDir.append("/") ;
3710  itosConv << msg->grpIndx ;
3711  baseDir += itosConv.str() ;
3712 
3713  if (stat(msg->baseDir, &info) != 0 ) {
3714  CkPrintf( "Node %d cannot access directory %s\n",
3715  CkMyPe(), baseDir.c_str() );
3716  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
3717  }
3718  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
3719  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
3720  int retVal = mkdir(baseDir.c_str(), S_IRWXU);
3721  }
3722 
3723  #ifdef DEBUG_QM
3724  Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
3725  #endif
3726 
3727  inputFileName.clear();
3728  inputFileName.append(baseDir.c_str()) ;
3729  inputFileName.append("/qmmm_") ;
3730  inputFileName += itosConv.str() ;
3731  inputFileName.append(".input") ;
3732 
3733  // Opens file for coordinate and parameter input
3734  inputFile = fopen(inputFileName.c_str(),"w");
3735  if ( ! inputFile ) {
3736  iout << iERROR << "Could not open input file for writing: "
3737  << inputFileName << "\n" << endi ;
3738  NAMD_err("Error writing QM input file.");
3739  }
3740 
3741  // Builds the command that will be executed
3742  qmCommand.clear();
3743  qmCommand.append("cd ");
3744  qmCommand.append(baseDir);
3745  qmCommand.append(" ; ");
3746  qmCommand.append(msg->execPath) ;
3747  qmCommand.append(" ") ;
3748  qmCommand.append(inputFileName) ;
3749 
3750  // Build the redirect file name we need for the screen output.
3751  // That's the only place where we can get partial charges for QM atoms.
3752  tmpRedirectFileName = inputFileName ;
3753  tmpRedirectFileName.append(".TmpOut") ;
3754 
3755  qmCommand.append(" > ") ;
3756  qmCommand.append(tmpRedirectFileName) ;
3757 
3758  // Builds the file name where orca will place the gradient
3759  // This will be relative to the input file
3760  outputFileName = inputFileName ;
3761  outputFileName.append(".engrad") ;
3762 
3763  // Builds the file name where orca will place gradients acting on
3764  // point charges
3765  pcGradFilename = inputFileName ;
3766  pcGradFilename.append(".pcgrad") ;
3767 
3768  if (msg->numAllPntChrgs) {
3769  // Builds the file name where we will write the point charges.
3770  pntChrgFileName = inputFileName ;
3771  pntChrgFileName.append(".pntchrg") ;
3772 
3773  pcGradFileName = inputFileName;
3774  pcGradFileName.append(".pcgrad");
3775 
3776  chrgFile = fopen(pntChrgFileName.c_str(),"w");
3777  if ( ! chrgFile ) {
3778  iout << iERROR << "Could not open charge file for writing: "
3779  << pntChrgFileName << "\n" << endi ;
3780  NAMD_err("Error writing point charge file.");
3781  }
3782 
3783  int numPntChrgs = 0;
3784  for (int i=0; i<msg->numAllPntChrgs; i++ ) {
3785  if (pcP[i].type != QMPCTYPE_IGNORE)
3786  numPntChrgs++;
3787  }
3788 
3789  iret = fprintf(chrgFile,"%d\n", numPntChrgs);
3790  if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
3791  }
3792 
3793  // writes configuration lines to the input file.
3794  std::stringstream ss(msg->configLines) ;
3795  std::string confLineStr;
3796  while (std::getline(ss, confLineStr) ) {
3797  confLineStr.append("\n");
3798  iret = fprintf(inputFile,confLineStr.c_str());
3799  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3800  }
3801 
3802  if (msg->numAllPntChrgs) {
3803  iret = fprintf(inputFile,"%%pointcharges \"%s\"\n", pntChrgFileName.c_str());
3804  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3805  }
3806 
3807  iret = fprintf(inputFile,"\n\n%%coords\n CTyp xyz\n");
3808  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3809 
3810  iret = fprintf(inputFile," Charge %f\n",msg->charge);
3811  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3812 
3813  iret = fprintf(inputFile," Mult %f\n",msg->multiplicity);
3814  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3815 
3816  iret = fprintf(inputFile," Units Angs\n coords\n\n");
3817  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3818 
3819  DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " <<
3820  inputFileName.c_str() << " and " << msg->numAllPntChrgs <<
3821  " point charges in file " << pntChrgFileName.c_str() << "\n");
3822 
3823  // write QM and dummy atom coordinates to input file.
3824  QMAtomData *atmP = msg->data ;
3825  for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
3826 
3827  double x = atmP->position.x;
3828  double y = atmP->position.y;
3829  double z = atmP->position.z;
3830 
3831  iret = fprintf(inputFile," %s %f %f %f\n",
3832  atmP->element,x,y,z);
3833  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3834 
3835  }
3836 
3837  iret = fprintf(inputFile," end\nend\n");
3838  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3839 
3840  if (msg->numAllPntChrgs) {
3841  // Write point charges to file.
3842  pcP = msg->data + msg->numAllAtoms ;
3843  for ( size_t j=0; j < msg->numAllPntChrgs; j++, ++pcP) {
3844 
3845  if (pcP->type == QMPCTYPE_IGNORE)
3846  continue;
3847 
3848  double charge = pcP->charge;
3849 
3850  double x = pcP->position.x;
3851  double y = pcP->position.y;
3852  double z = pcP->position.z;
3853 
3854  iret = fprintf(chrgFile,"%f %f %f %f\n",
3855  charge,x,y,z);
3856  if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
3857  }
3858 
3859  fclose(chrgFile);
3860  }
3861 
3862  DebugM(4,"Closing input and charge file\n");
3863  fclose(inputFile);
3864 
3865  if (msg->prepProcOn) {
3866 
3867  std::string prepProc(msg->prepProc) ;
3868  prepProc.append(" ") ;
3869  prepProc.append(inputFileName) ;
3870  iret = system(prepProc.c_str());
3871  if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
3872  if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
3873  }
3874 
3875  // runs QM command
3876  DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
3877  iret = system(qmCommand.c_str());
3878 
3879  if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
3880  if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
3881 
3882  if (msg->secProcOn) {
3883 
3884  std::string secProc(msg->secProc) ;
3885  secProc.append(" ") ;
3886  secProc.append(inputFileName) ;
3887  itosConv.str("");
3888  itosConv.clear() ;
3889  itosConv << msg->timestep ;
3890  secProc.append(" ") ;
3891  secProc += itosConv.str() ;
3892 
3893  iret = system(secProc.c_str());
3894  if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
3895  if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
3896  }
3897 
3898  // remove coordinate file
3899 // iret = remove(inputFileName);
3900 // if ( iret ) { NAMD_err(0); }
3901 
3902  // remove coordinate file
3903 // iret = remove(pntChrgFileName);
3904 // if ( iret ) { NAMD_err(0); }
3905 
3906  // opens output file
3907  DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
3908  outputFile = fopen(outputFileName.c_str(),"r");
3909  if ( ! outputFile ) {
3910  iout << iERROR << "Could not find QM output file!\n" << endi;
3911  NAMD_err(0);
3912  }
3913 
3914  // Resets the pointers.
3915  atmP = msg->data ;
3916  pcP = msg->data + msg->numAllAtoms ;
3917 
3918  // Allocates the return message.
3919  QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
3920  resMsg->grpIndx = msg->grpIndx;
3921  resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
3922  resMsg->energyOrig = 0;
3923  resMsg->energyCorr = 0;
3924  for ( int k=0; k<3; ++k )
3925  for ( int l=0; l<3; ++l )
3926  resMsg->virial[k][l] = 0;
3927  QMForce *resForce = resMsg->force ;
3928 
3929  // We split the data into two pointers so we can "skip" the dummy atoms
3930  // which have no representation in NAMD.
3931  for (int i=0; i<resMsg->numForces; i++) {
3932  resForce[i].force = 0;
3933  resForce[i].charge = 0 ;
3934  if (i < msg->numQMAtoms) {
3935  // We use the replace field to indicate QM atoms,
3936  // which will have charge information.
3937  resForce[i].replace = 1 ;
3938  resForce[i].id = atmP->id;
3939  atmP++;
3940  }
3941  else {
3942  // We use the replace field to indicate QM atoms,
3943  // which will have charge information.
3944  resForce[i].replace = 0 ;
3945  resForce[i].id = pcP->id;
3946  pcP++;
3947  }
3948  }
3949 
3950  // Resets the pointers.
3951  atmP = msg->data ;
3952  pcP = msg->data + msg->numAllAtoms ;
3953 
3954  size_t atmIndx = 0, gradCount = 0;
3955  int numAtomsInFile = 0;
3956 
3957  // Reads the data form the output file created by the QM software.
3958  // Gradients over the QM atoms, and Charges for QM atoms will be read.
3959 
3960  // skip lines before number of atoms
3961  for (int i = 0; i < 3; i++) {
3962  fgets(line, lineLen, outputFile);
3963  }
3964 
3965  iret = fscanf(outputFile,"%d\n", &numAtomsInFile);
3966  if ( iret != 1 ) {
3967  NAMD_die("Error reading QM forces file.");
3968  }
3969 
3970  #ifdef DEBUG_QM
3971  if(numAtomsInFile != msg->numAllAtoms) {
3972  NAMD_die("Error reading QM forces file. Number of atoms in QM output\
3973  does not match the expected.");
3974  }
3975  #endif
3976 
3977  // skip lines before energy
3978  for (int i = 0; i < 3; i++) {
3979  fgets(line, lineLen, outputFile);
3980  }
3981 
3982  iret = fscanf(outputFile,"%lf\n", &resMsg->energyOrig);
3983  if ( iret != 1 ) {
3984  NAMD_die("Error reading QM forces file.");
3985  }
3986 // iout << "Energy in step (Hartree): " << resMsg->energyOrig << "\n" << endi;
3987  // All energies are given in Eh (Hartree)
3988  // NAMD needs energies in kcal/mol
3989  // The conversion factor is 627.509469
3990  resMsg->energyOrig *= 627.509469;
3991 // iout << "Energy in step (Kcal/mol): " << resMsg->energyOrig << "\n" << endi;
3992 
3993  resMsg->energyCorr = resMsg->energyOrig;
3994 
3995  // skip lines before gradient
3996  for (int i = 0; i < 3; i++) {
3997  fgets(line, lineLen, outputFile) ;
3998  }
3999 
4000  // Break the loop when we find all gradients for QM atoms,
4001  // and keep the next gradient lines from being read, if there
4002  // are more. Following gradients, if present, will refer to link
4003  // atoms.
4004  atmIndx = 0;
4005  gradCount = 0;
4006  for ( size_t i=0; i<msg->numAllAtoms*3; ++i ) {
4007 
4008  iret = fscanf(outputFile,"%lf\n", &gradient[gradCount]);
4009  if ( iret != 1 ) { NAMD_die("Error reading QM forces file."); }
4010 
4011  if (gradCount == 2){
4012 
4013  // All gradients are given in Eh/a0 (Hartree over Bohr radius)
4014  // NAMD needs forces in kcal/mol/angstrons
4015  // The conversion factor is 627.509469/0.529177 = 1185.82151
4016  if (atmIndx < msg->numQMAtoms ) {
4017  resForce[atmIndx].force.x = -1*gradient[0]*1185.82151;
4018  resForce[atmIndx].force.y = -1*gradient[1]*1185.82151;
4019  resForce[atmIndx].force.z = -1*gradient[2]*1185.82151;
4020  }
4021 
4022  // If we are reading forces applied on Dummy atoms,
4023  // place them on the appropriate QM and MM atoms to conserve energy.
4024 
4025  // This implementation was based on the description in
4026  // DOI: 10.1002/jcc.20857 : Equations 30 to 32
4027  if ( msg->numQMAtoms <= atmIndx &&
4028  atmIndx < msg->numAllAtoms ) {
4029  // The dummy atom points to the QM atom to which it is bound.
4030  // The QM atom and the MM atom (in a QM-MM bond) point to each other.
4031  int qmInd = atmP[atmIndx].bountToIndx ;
4032  int mmInd = atmP[qmInd].bountToIndx ;
4033 
4034  Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
4035  Real mmqmDist = dir.length() ;
4036 
4037  Real linkDist = Vector(atmP[atmIndx].position - atmP[qmInd].position).length() ;
4038 
4039  Force mmForce(0), qmForce(0), linkForce(gradient[0], gradient[1], gradient[2]);
4040  linkForce *= -1*1185.82151;
4041 
4042  Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
4043  // Unit vectors
4044  Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
4045  Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
4046  Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
4047  Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
4048 
4049  qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv +
4050  (xDelta)*base) )*xuv;
4051 
4052  qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv +
4053  (yDelta)*base) )*yuv;
4054 
4055  qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv +
4056  (zDelta)*base) )*zuv;
4057 
4058 
4059  mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
4060  (xDelta)*base) )*xuv;
4061 
4062  mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
4063  (yDelta)*base) )*yuv;
4064 
4065  mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
4066  (zDelta)*base) )*zuv;
4067 
4068  resForce[qmInd].force += qmForce;
4069  resForce[msg->numQMAtoms + mmInd].force += mmForce;
4070  }
4071 
4072  gradCount = 0;
4073  atmIndx++ ;
4074  } else
4075  gradCount++ ;
4076 
4077  }
4078 
4079  fclose(outputFile);
4080 
4081  // In case charges are not to be read form the QM software output,
4082  // we load the origianl atom charges.
4083  if (msg->qmAtmChrgMode == QMCHRGNONE) {
4084  int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
4085  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
4086  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
4087 
4088  atmIndx = 0 ;
4089  for (; atmIndx < msg->numQMAtoms; atmIndx++) {
4090 
4091  // Loops over all QM atoms (in all QM groups) comparing their global indices
4092  for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
4093 
4094  if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
4095 
4096  atmP[atmIndx].charge = qmAtmChrg[qmIter];
4097  resForce[atmIndx].charge = qmAtmChrg[qmIter];
4098 
4099  break;
4100  }
4101 
4102  }
4103 
4104  }
4105  }
4106  else {
4107  // opens redirected output file
4108  outputFile = fopen(tmpRedirectFileName.c_str(),"r");
4109  if ( ! outputFile ) {
4110  iout << iERROR << "Could not find Redirect output file:"
4111  << tmpRedirectFileName << std::endl << endi;
4112  NAMD_err(0);
4113  }
4114 
4115  DebugM(4,"Opened tmeporary output for charge reading: " << tmpRedirectFileName.c_str() << "\n");
4116 
4117  int lineState = 0;
4118  char result[20] ;
4119  while ( fgets(line, lineLen, outputFile) != NULL){
4120 
4121  // We loop over the lines of the output file untill we find
4122  // the line that initiates the charges lines. Then
4123  // we read all lines and expect to get one charge
4124  // per line, untill we find a line that sums all charges.
4125 
4126  switch (msg->qmAtmChrgMode) {
4127 
4128  case QMCHRGMULLIKEN:
4129  {
4130 
4131  if ( strstr(line,"MULLIKEN ATOMIC CHARGES") != NULL ) {
4132  lineState = 1;
4133  atmIndx = 0;
4134 
4135  // Now we expect the following line to have a series of dashes
4136  // and the folowing lines to have atom charges (among other info)
4137  continue;
4138  }
4139 
4140  if ( strstr(line,"Sum of atomic charges") != NULL ) {
4141 
4142  // Now that all charges have been read, grabs the total charge
4143  // just for fun... and checking purposes.
4144  #ifdef DEBUG_QM
4145  strncpy(result, line + 31,12) ;
4146  result[12] = '\0';
4147 
4148  DebugM(4,"Total charge of QM region calculated by ORCA is: "
4149  << atof(result) << std::endl )
4150  #endif
4151 
4152  // Nothing more to read from the output file
4153  break;
4154  }
4155 
4156  // Line state 1 means we have to skip ONLY one line.
4157  if (lineState == 1) {
4158  lineState = 2;
4159  continue;
4160  }
4161 
4162  // Line state 2 means we have to read the line and grab the info.
4163  if (lineState == 2) {
4164 
4165  strncpy(result, line+8,12) ;
4166  result[12] = '\0';
4167 
4168  Real localCharge = atof(result);
4169 
4170  // If we are reading charges from QM atoms, store them.
4171  if (atmIndx < msg->numQMAtoms ) {
4172  atmP[atmIndx].charge = localCharge;
4173  resForce[atmIndx].charge = localCharge;
4174  }
4175 
4176  // If we are reading charges from Dummy atoms,
4177  // place the on the appropriate QM atom.
4178  if ( msg->numQMAtoms <= atmIndx &&
4179  atmIndx < msg->numAllAtoms ) {
4180  int qmInd = atmP[atmIndx].bountToIndx ;
4181  atmP[qmInd].charge += localCharge;
4182  resForce[qmInd].charge += localCharge;
4183  }
4184 
4185  atmIndx++ ;
4186 
4187  // If we found all charges for QM atoms, change the lineState
4188  // untill we reach the "total charge" line.
4189  if (atmIndx == msg->numAllAtoms ) {
4190  lineState = 0;
4191  }
4192 
4193  continue;
4194  }
4195 
4196  } break ;
4197 
4198  case QMCHRGCHELPG :
4199  {
4200 
4201  if ( strstr(line,"CHELPG Charges") != NULL ) {
4202  lineState = 1;
4203  atmIndx = 0;
4204 
4205  // Now we expect the following line to have a series of dashes
4206  // and the folowing lines to have atom charges (among other info)
4207  continue;
4208  }
4209 
4210  if ( strstr(line,"Total charge") != NULL ) {
4211 
4212  // Now that all charges have been read, grabs the total charge
4213  // just for fun... and checking purposes.
4214  #ifdef DEBUG_QM
4215  strncpy(result, line + 14,13) ;
4216  result[13] = '\0';
4217 
4218  DebugM(4,"Total charge of QM region calculated by ORCA is: "
4219  << atof(result) << std::endl )
4220  #endif
4221 
4222  // Nothing more to read from the output file
4223  break;
4224  }
4225 
4226  // Line state 1 means we have to skip ONLY one line.
4227  if (lineState == 1) {
4228  lineState = 2;
4229  continue;
4230  }
4231 
4232  // Line state 2 means we have to read the line and grab the info.
4233  if (lineState == 2) {
4234 
4235  strncpy(result, line+12,15) ;
4236  result[15] = '\0';
4237 
4238  Real localCharge = atof(result);
4239 
4240  // If we are reading charges from QM atoms, store them.
4241  if (atmIndx < msg->numQMAtoms ) {
4242  atmP[atmIndx].charge = localCharge;
4243  resForce[atmIndx].charge = localCharge;
4244  }
4245 
4246  // If we are reading charges from Dummy atoms,
4247  // place the on the appropriate QM atom.
4248  if ( msg->numQMAtoms <= atmIndx &&
4249  atmIndx < msg->numAllAtoms ) {
4250  int qmInd = atmP[atmIndx].bountToIndx ;
4251  atmP[qmInd].charge += localCharge;
4252  resForce[qmInd].charge += localCharge;
4253  }
4254 
4255  atmIndx++ ;
4256 
4257  // If we found all charges for QM atoms, we ignore the following line
4258  // untill we reach the "total charge" line.
4259  if (atmIndx == msg->numAllAtoms ) {
4260  lineState = 1;
4261  }
4262 
4263  continue;
4264  }
4265 
4266  } break;
4267  }
4268  }
4269 
4270  fclose(outputFile);
4271  }
4272 
4273  // remove force file
4274 // DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
4275 // iret = remove(outputFileName);
4276 // if ( iret ) { NAMD_err(0); }
4277 
4278 
4279  DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
4280 
4281  atmP = msg->data ;
4282  pcP = msg->data + msg->numAllAtoms ;
4283 
4284  // The initial point charge index for the force message is the number of QM
4285  // atoms, since the dummy atoms have no representation in NAMD
4286  int pcIndx = msg->numQMAtoms;
4287 
4288  // If there are no point charges, orca will not create a "pcgrad" file.
4289  if (msg->numAllPntChrgs > 0) {
4290 
4291  // opens redirected output file
4292  outputFile = fopen(pcGradFileName.c_str(),"r");
4293  if ( ! outputFile ) {
4294  iout << iERROR << "Could not find PC gradient output file:"
4295  << pcGradFileName << std::endl << endi;
4296  NAMD_err(0);
4297  }
4298 
4299  DebugM(4,"Opened pc output for gradient reading: " << pcGradFileName.c_str() << "\n");
4300 
4301  int pcgradNumPC = 0, readPC = 0;
4302  iret = fscanf(outputFile,"%d\n", &pcgradNumPC );
4303  if ( iret != 1 ) {
4304  NAMD_die("Error reading PC forces file.");
4305  }
4306  DebugM(4,"Found in pc gradient output: " << pcgradNumPC << " point charge grads.\n");
4307 
4308  // We loop over point charges, reading the total electrostatic force
4309  // applied on them by the QM region.
4310  // We redistribute the forces applied over virtual point
4311  // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
4312  for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
4313 
4314  Force totalForce(0);
4315 
4316  // No force was applied to the QM region due to this charge, since it
4317  // was skipped when writing down point charges to the QM software, so it
4318  // does not receive any back from the QM region. It must be an MM1 atom
4319  // from a QM-MM bond.
4320  if (pcP[i].type == QMPCTYPE_IGNORE)
4321  continue;
4322 
4323  fgets(line, lineLen, outputFile);
4324 
4325  iret = sscanf(line,"%lf %lf %lf\n", &totalForce[0], &totalForce[1], &totalForce[2] );
4326  if ( iret != 3 ) {
4327  NAMD_die("Error reading PC forces file.");
4328  }
4329  // All gradients are given in Eh/a0 (Hartree over Bohr radius)
4330  // NAMD needs forces in kcal/mol/angstrons
4331  // The conversion factor is 627.509469/0.529177 = 1185.82151
4332  totalForce *= -1185.82151;
4333 
4334  if (pcP[i].type == QMPCTYPE_CLASSICAL) {
4335  // Checking pcP was not a QM atom in another region
4336  // If so the interaction is accounted in the other region
4337  if (qmAtomGroup[pcP[i].id] == 0) {
4338  // If we already ignored MM1 charges, we take all other
4339  // non-virtual charges and apply forces directly to them.
4340  resForce[pcIndx].force += totalForce;
4341  }
4342  }
4343  else {
4344  // If we are handling virtual PC, we distribute the force over
4345  // MM1 and MM2.
4346 
4347  Force mm1Force(0), mm2Force(0);
4348 
4349  // Virtual PC are bound to MM2.
4350  int mm2Indx = pcP[i].bountToIndx;
4351  // MM2 charges are bound to MM1.
4352  int mm1Indx = pcP[mm2Indx].bountToIndx;
4353 
4354  Real Cq = pcP[i].dist;
4355 
4356  mm1Force = (1-Cq)*totalForce ;
4357  mm2Force = Cq*totalForce ;
4358 
4359  resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
4360  resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
4361  }
4362 
4363  }
4364 
4365  fclose(outputFile);
4366  }
4367 
4368  delete [] line;
4369 
4370  // Adjusts forces from PME, canceling contributions from the QM and
4371  // direct Coulomb forces calculated here.
4372  if (msg->PMEOn) {
4373 
4374  DebugM(1,"Correcting forces and energy for PME.\n");
4375 
4376  ewaldcof = msg->PMEEwaldCoefficient;
4377  BigReal TwoBySqrtPi = 1.12837916709551;
4378  pi_ewaldcof = TwoBySqrtPi * ewaldcof;
4379 
4380  for (size_t i=0; i < msg->numQMAtoms; i++) {
4381 
4382  BigReal p_i_charge = atmP[i].charge ;
4383  Position pos_i = atmP[i].position ;
4384 
4385  for (size_t j=i+1; j < msg->numQMAtoms; j++) {
4386 
4387  BigReal p_j_charge = atmP[j].charge ;
4388 
4389  Position pos_j = atmP[j].position ;
4390 
4391  BigReal r = Vector(pos_i - pos_j).length() ;
4392 
4393  BigReal tmp_a = r * ewaldcof;
4394  BigReal tmp_b = erfc(tmp_a);
4395  BigReal corr_energy = tmp_b;
4396  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
4397 
4398 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
4399  BigReal recip_energy = (1-tmp_b)/r;
4400 
4401 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
4402  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
4403 
4404  const BigReal kq_i = p_i_charge * constants;
4405 
4406  // Final force and energy correction for this pair of atoms.
4407  BigReal energy = kq_i * p_j_charge * recip_energy ;
4408 
4409  Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
4410 
4411  // The force is *subtracted* from the total force acting on
4412  // both atoms. The sign on fixForce corrects the orientation
4413  // of the subtracted force.
4414 // DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
4415 // << std::endl);
4416 // DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
4417 // << std::endl);
4418 // DebugM(4,"Force correction: " << fixForce << std::endl);
4419  resForce[i].force -= fixForce ;
4420  resForce[j].force -= -1*fixForce;
4421 
4422  // The energy is *subtracted* from the total energy calculated here.
4423  resMsg->energyCorr -= energy;
4424  }
4425  }
4426 
4427  pcIndx = msg->numQMAtoms;
4428  // We only loop over point charges from real atoms, ignoring the ones
4429  // created to handle QM-MM bonds.
4430  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
4431 
4432  BigReal p_i_charge = pcPpme[i].charge;
4433  Position pos_i = pcPpme[i].position ;
4434 
4435  const BigReal kq_i = p_i_charge * constants;
4436 
4437  Force fixForce = 0;
4438 
4439  for (size_t j=0; j<msg->numQMAtoms; ++j ) {
4440 
4441  BigReal p_j_charge = atmP[j].charge ;
4442 
4443  Position pos_j = atmP[j].position ;
4444 
4445  BigReal r = Vector(pos_i - pos_j).length() ;
4446 
4447  BigReal tmp_a = r * ewaldcof;
4448  BigReal tmp_b = erfc(tmp_a);
4449  BigReal corr_energy = tmp_b;
4450  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
4451 
4452 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
4453  BigReal recip_energy = (1-tmp_b)/r;
4454 
4455 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
4456  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
4457 
4458  // Final force and energy correction for this pair of atoms.
4459  BigReal energy = kq_i * p_j_charge * recip_energy ;
4460 
4461  fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
4462 
4463  // The energy is *subtracted* from the total energy calculated here.
4464  resMsg->energyCorr -= energy;
4465 
4466  }
4467 
4468  // The force is *subtracted* from the total force acting on
4469  // the point charge.
4470 // DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
4471 // << std::endl);
4472 // DebugM(4,"Force correction: " << fixForce << std::endl);
4473  resForce[pcIndx].force -= kq_i*fixForce ;
4474  }
4475 
4476  }
4477 
4478  DebugM(1,"Determining virial...\n");
4479 
4480  // Calculates the virial contribution form the forces on QM atoms and
4481  // point charges calculated here.
4482  for (size_t i=0; i < msg->numQMAtoms; i++) {
4483  // virial used by NAMD is -'ve of normal convention, so reverse it!
4484  // virial[i][j] in file should be sum of -1 * f_i * r_j
4485  for ( int k=0; k<3; ++k )
4486  for ( int l=0; l<3; ++l )
4487  resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
4488  }
4489 
4490  pcIndx = msg->numQMAtoms; // Index in the real PC force array.
4491  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
4492  // virial used by NAMD is -'ve of normal convention, so reverse it!
4493  // virial[i][j] in file should be sum of -1 * f_i * r_j
4494  for ( int k=0; k<3; ++k )
4495  for ( int l=0; l<3; ++l )
4496  resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
4497  }
4498 
4499  DebugM(1,"End of QM processing. Sending result message.\n");
4500 
4501  // Send message to rank zero with results.
4502  QMProxy[0].recvQMRes(resMsg);
4503 
4504  if (msg->switching && msg->PMEOn)
4505  delete [] pcPpme;
4506 
4507  delete msg;
4508  return ;
4509 }
4510 
4512  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
4513  FILE *inputFile,*outputFile;
4514  int iret;
4515 
4516  std::string qmCommand, inputFileName, outputFileName;
4517 
4518  // For coulomb calculation
4519  BigReal constants = msg->constants;
4520 
4521  DebugM(4,"Running USER DEFINED SOFTWARE on PE " << CkMyPe() << std::endl);
4522 
4523  QMAtomData *pcP = msg->data + msg->numAllAtoms ;
4524 
4525  // We create a pointer for PME correction, which may point to
4526  // a copy of the original point charge array, with unchanged charges,
4527  // or points to the original array in case no switching or charge
4528  // scheme is being used.
4529  QMAtomData *pcPpme = NULL;
4530  if (msg->switching) {
4531 
4532  if (msg->PMEOn)
4533  pcPpme = new QMAtomData[msg->numRealPntChrgs];
4534 
4535  pntChrgSwitching(msg, pcPpme) ;
4536  } else {
4537  pcPpme = pcP;
4538  }
4539 
4540  int retVal = 0;
4541  struct stat info;
4542 
4543  // For each QM group, create a subdirectory where files will be palced.
4544  std::string baseDir(msg->baseDir);
4545  std::ostringstream itosConv ;
4546  if ( CmiNumPartitions() > 1 ) {
4547  baseDir.append("/") ;
4548  itosConv << CmiMyPartition() ;
4549  baseDir += itosConv.str() ;
4550  itosConv.str("");
4551  itosConv.clear() ;
4552 
4553  if (stat(msg->baseDir, &info) != 0 ) {
4554  CkPrintf( "Node %d cannot access directory %s\n",
4555  CkMyPe(), baseDir.c_str() );
4556  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
4557  }
4558  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
4559  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
4560  retVal = mkdir(baseDir.c_str(), S_IRWXU);
4561  }
4562 
4563  }
4564  baseDir.append("/") ;
4565  itosConv << msg->grpIndx ;
4566  baseDir += itosConv.str() ;
4567 
4568  if (stat(msg->baseDir, &info) != 0 ) {
4569  CkPrintf( "Node %d cannot access directory %s\n",
4570  CkMyPe(), baseDir.c_str() );
4571  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
4572  }
4573  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
4574  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
4575  int retVal = mkdir(baseDir.c_str(), S_IRWXU);
4576  }
4577 
4578  #ifdef DEBUG_QM
4579  Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
4580  #endif
4581 
4582  inputFileName.clear();
4583  inputFileName.append(baseDir.c_str()) ;
4584  inputFileName.append("/qmmm_") ;
4585  inputFileName += itosConv.str() ;
4586  inputFileName.append(".input") ;
4587 
4588  // Opens file for coordinate and parameter input
4589  inputFile = fopen(inputFileName.c_str(),"w");
4590  if ( ! inputFile ) {
4591  iout << iERROR << "Could not open input file for writing: "
4592  << inputFileName << "\n" << endi ;
4593  NAMD_err("Error writing QM input file.");
4594  }
4595 
4596  // Builds the command that will be executed
4597  qmCommand.clear();
4598  qmCommand.append("cd ");
4599  qmCommand.append(baseDir);
4600  qmCommand.append(" ; ");
4601  qmCommand.append(msg->execPath) ;
4602  qmCommand.append(" ") ;
4603  qmCommand.append(inputFileName) ;
4604 
4605  // Builds the file name where orca will place the gradient
4606  // This will be relative to the input file
4607  outputFileName = inputFileName ;
4608  outputFileName.append(".result") ;
4609 
4610  int numPntChrgs = 0;
4611  for (int i=0; i<msg->numAllPntChrgs; i++ ) {
4612  if (pcP[i].type != QMPCTYPE_IGNORE)
4613  numPntChrgs++;
4614  }
4615 
4616  iret = fprintf(inputFile,"%d %d\n",msg->numAllAtoms, numPntChrgs);
4617  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
4618 
4619  DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " <<
4620  inputFileName.c_str() << " and " << msg->numAllPntChrgs <<
4621  " point charges." << std::endl);
4622 
4623  // write QM and dummy atom coordinates to input file.
4624  QMAtomData *atmP = msg->data ;
4625  for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
4626 
4627  double x = atmP->position.x;
4628  double y = atmP->position.y;
4629  double z = atmP->position.z;
4630 
4631  iret = fprintf(inputFile,"%f %f %f %s\n",
4632  x,y,z,atmP->element);
4633  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
4634 
4635  }
4636 
4637  int numWritenPCs = 0;
4638  // Write point charges to file.
4639  pcP = msg->data + msg->numAllAtoms ;
4640  for ( size_t j=0; j < msg->numAllPntChrgs; j++, ++pcP) {
4641 
4642  if (pcP->type == QMPCTYPE_IGNORE)
4643  continue;
4644 
4645  double charge = pcP->charge;
4646 
4647  double x = pcP->position.x;
4648  double y = pcP->position.y;
4649  double z = pcP->position.z;
4650 
4651  iret = fprintf(inputFile,"%f %f %f %f\n",
4652  x,y,z,charge);
4653  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
4654 
4655  numWritenPCs++;
4656  }
4657 
4658  DebugM(4,"Closing input file\n");
4659  fclose(inputFile);
4660 
4661  if (msg->prepProcOn) {
4662 
4663  std::string prepProc(msg->prepProc) ;
4664  prepProc.append(" ") ;
4665  prepProc.append(inputFileName) ;
4666  iret = system(prepProc.c_str());
4667  if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
4668  if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
4669  }
4670 
4671  // runs QM command
4672  DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
4673  iret = system(qmCommand.c_str());
4674 
4675  if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
4676  if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
4677 
4678  if (msg->secProcOn) {
4679 
4680  std::string secProc(msg->secProc) ;
4681  secProc.append(" ") ;
4682  secProc.append(inputFileName) ;
4683  itosConv.str("");
4684  itosConv.clear() ;
4685  itosConv << msg->timestep ;
4686  secProc.append(" ") ;
4687  secProc += itosConv.str() ;
4688 
4689  iret = system(secProc.c_str());
4690  if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
4691  if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
4692  }
4693 
4694  // remove coordinate file
4695 // iret = remove(inputFileName);
4696 // if ( iret ) { NAMD_err(0); }
4697 
4698  // remove coordinate file
4699 // iret = remove(pntChrgFileName);
4700 // if ( iret ) { NAMD_err(0); }
4701 
4702  // opens output file
4703  DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
4704  outputFile = fopen(outputFileName.c_str(),"r");
4705  if ( ! outputFile ) {
4706  iout << iERROR << "Could not find QM output file!\n" << endi;
4707  NAMD_err(0);
4708  }
4709 
4710  // Resets the pointers.
4711  atmP = msg->data ;
4712  pcP = msg->data + msg->numAllAtoms ;
4713 
4714  // Allocates the return message.
4715  QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
4716  resMsg->grpIndx = msg->grpIndx;
4717  resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
4718  resMsg->energyOrig = 0;
4719  resMsg->energyCorr = 0;
4720  for ( int k=0; k<3; ++k )
4721  for ( int l=0; l<3; ++l )
4722  resMsg->virial[k][l] = 0;
4723  QMForce *resForce = resMsg->force ;
4724 
4725  // We split the data into two pointers so we can "skip" the dummy atoms
4726  // which have no representation in NAMD.
4727  for (int i=0; i<resMsg->numForces; i++) {
4728  resForce[i].force = 0;
4729  resForce[i].charge = 0 ;
4730  if (i < msg->numQMAtoms) {
4731  // We use the replace field to indicate QM atoms,
4732  // which will have charge information.
4733  resForce[i].replace = 1 ;
4734  resForce[i].id = atmP->id;
4735  atmP++;
4736  }
4737  else {
4738  // We use the replace field to indicate QM atoms,
4739  // which will have charge information.
4740  resForce[i].replace = 0 ;
4741  resForce[i].id = pcP->id;
4742  pcP++;
4743  }
4744  }
4745 
4746  // Resets the pointers.
4747  atmP = msg->data ;
4748  pcP = msg->data + msg->numAllAtoms ;
4749 
4750  // Reads the data form the output file created by the QM software.
4751  // Gradients over the QM atoms, and Charges for QM atoms will be read.
4752 
4753  // Number of point charges for which we will receive forces.
4754  int usrPCnum = 0;
4755  const size_t lineLen = 256;
4756  char *line = new char[lineLen];
4757 
4758  fgets(line, lineLen, outputFile);
4759 
4760 // iret = fscanf(outputFile,"%lf %d\n", &resMsg->energyOrig, &usrPCnum);
4761  iret = sscanf(line,"%lf %i\n", &resMsg->energyOrig, &usrPCnum);
4762  if ( iret < 1 ) {
4763  NAMD_die("Error reading energy from QM results file.");
4764  }
4765 
4766  resMsg->energyCorr = resMsg->energyOrig;
4767 
4768  if (iret == 2 && numWritenPCs != usrPCnum) {
4769  iout << iERROR << "Number of point charges does not match what was provided!\n" << endi ;
4770  NAMD_die("Error reading QM results file.");
4771  }
4772 
4773  size_t atmIndx;
4774  double localForce[3];
4775  double localCharge;
4776  for (atmIndx = 0; atmIndx < msg->numAllAtoms; atmIndx++) {
4777 
4778  iret = fscanf(outputFile,"%lf %lf %lf %lf\n",
4779  localForce+0,
4780  localForce+1,
4781  localForce+2,
4782  &localCharge);
4783  if ( iret != 4 ) {
4784  NAMD_die("Error reading forces and charge from QM results file.");
4785  }
4786 
4787  // If we are reading charges and forces on QM atoms, store
4788  // them directly.
4789  if (atmIndx < msg->numQMAtoms ) {
4790 
4791  resForce[atmIndx].force.x = localForce[0];
4792  resForce[atmIndx].force.y = localForce[1];
4793  resForce[atmIndx].force.z = localForce[2];
4794 
4795  atmP[atmIndx].charge = localCharge;
4796  resForce[atmIndx].charge = localCharge;
4797  }
4798 
4799  // If we are reading charges and forces applied on Dummy atoms,
4800  // place them on the appropriate QM and MM atoms to conserve energy.
4801 
4802  // This implementation was based on the description in
4803  // DOI: 10.1002/jcc.20857 : Equations 30 to 32
4804  if ( msg->numQMAtoms <= atmIndx &&
4805  atmIndx < msg->numAllAtoms ) {
4806 
4807  // We keep the calculated charge in the dummy atom
4808  // so that we can calculate electrostatic forces later on.
4809  atmP[atmIndx].charge = localCharge;
4810 
4811  // If we are reading charges from Dummy atoms,
4812  // place them on the appropriate QM atom.
4813  // The dummy atom points to the QM atom to which it is bound.
4814  int qmInd = atmP[atmIndx].bountToIndx ;
4815  resForce[qmInd].charge += localCharge;
4816 
4817  // The dummy atom points to the QM atom to which it is bound.
4818  // The QM atom and the MM atom (in a QM-MM bond) point to each other.
4819  int mmInd = atmP[qmInd].bountToIndx ;
4820 
4821  Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
4822  Real mmqmDist = dir.length() ;
4823 
4824  Real linkDist = Vector(atmP[atmIndx].position -
4825  atmP[qmInd].position).length() ;
4826 
4827  Force mmForce(0), qmForce(0),
4828  linkForce(localForce[0], localForce[1], localForce[2]);
4829 
4830  Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
4831  // Unit vectors
4832  Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
4833  Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
4834  Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
4835  Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
4836 
4837  qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv +
4838  (xDelta)*base) )*xuv;
4839 
4840  qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv +
4841  (yDelta)*base) )*yuv;
4842 
4843  qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv +
4844  (zDelta)*base) )*zuv;
4845 
4846 
4847  mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
4848  (xDelta)*base) )*xuv;
4849 
4850  mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
4851  (yDelta)*base) )*yuv;
4852 
4853  mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
4854  (zDelta)*base) )*zuv;
4855 
4856  resForce[qmInd].force += qmForce;
4857  resForce[msg->numQMAtoms + mmInd].force += mmForce;
4858  }
4859  }
4860 
4861  // The initial point charge index for the force message is the number of QM
4862  // atoms, since the dummy atoms have no representation in NAMD
4863  int pcIndx = msg->numQMAtoms;
4864 
4865  if (usrPCnum > 0) {
4866  // We loop over point charges, reading the total electrostatic force
4867  // applied on them by the QM region.
4868  // We redistribute the forces applied over virtual point
4869  // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
4870  for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
4871 
4872  Force totalForce(0);
4873 
4874  // No force was applied to the QM region due to this charge, since it
4875  // was skipped when writing down point charges to the QM software, so it
4876  // does not receive any back from the QM region. It must be an MM1 atom
4877  // from a QM-MM bond.
4878  if (pcP[i].type == QMPCTYPE_IGNORE)
4879  continue;
4880 
4881  iret = fscanf(outputFile,"%lf %lf %lf\n",
4882  &totalForce[0], &totalForce[1], &totalForce[2]);
4883  if ( iret != 3 ) {
4884  NAMD_die("Error reading PC forces from QM results file.");
4885  }
4886 
4887  if (pcP[i].type == QMPCTYPE_CLASSICAL) {
4888  // Checking pcP was not a QM atom in another region
4889  // If so the interaction is accounted in the other region
4890  if (qmAtomGroup[pcP[i].id] == 0) {
4891  // If we already ignored MM1 charges, we take all other
4892  // non-virtual charges and apply forces directly to them.
4893  resForce[pcIndx].force += totalForce;
4894  }
4895  }
4896  else {
4897  // If we are handling virtual PC, we distribute the force over
4898  // MM1 and MM2.
4899 
4900  Force mm1Force(0), mm2Force(0);
4901 
4902  // Virtual PC are bound to MM2.
4903  int mm2Indx = pcP[i].bountToIndx;
4904  // MM2 charges are bound to MM1.
4905  int mm1Indx = pcP[mm2Indx].bountToIndx;
4906 
4907  Real Cq = pcP[i].dist;
4908 
4909  mm1Force = (1-Cq)*totalForce ;
4910  mm2Force = Cq*totalForce ;
4911 
4912  resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
4913  resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
4914  }
4915 
4916 
4917  }
4918  }
4919 
4920  fclose(outputFile);
4921  delete [] line;
4922 
4923  // In case charges are not to be read form the QM software output,
4924  // we load the origianl atom charges.
4925  if (msg->qmAtmChrgMode == QMCHRGNONE) {
4926  int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
4927  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
4928  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
4929 
4930  atmIndx = 0 ;
4931  for (; atmIndx < msg->numQMAtoms; atmIndx++) {
4932 
4933  // Loops over all QM atoms (in all QM groups) comparing their global indices
4934  for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
4935 
4936  if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
4937 
4938  atmP[atmIndx].charge = qmAtmChrg[qmIter];
4939  resForce[atmIndx].charge = qmAtmChrg[qmIter];
4940 
4941  break;
4942  }
4943 
4944  }
4945 
4946  }
4947  for (size_t j=msg->numQMAtoms; j<msg->numAllAtoms; ++j ) {
4948  atmP[j].charge = 0;
4949  }
4950  }
4951 
4952  // remove force file
4953 // DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
4954 // iret = remove(outputFileName);
4955 // if ( iret ) { NAMD_err(0); }
4956 
4957  if (usrPCnum == 0) {
4958  DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
4959 
4960  atmP = msg->data ;
4961  pcP = msg->data + msg->numAllAtoms ;
4962 
4963  // We only loop over point charges from real atoms, ignoring the ones
4964  // created to handle QM-MM bonds.
4965  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
4966 
4967  // No force was applied to the QM region due to this charge, so it
4968  // does not receive any back from the QM region. It must be an MM1 atom
4969  // from a QM-MM bond.
4970  if (pcP[i].type == QMPCTYPE_IGNORE)
4971  continue;
4972 
4973  Force totalForce(0);
4974 
4975  BigReal pntCharge = pcP[i].charge;
4976 
4977  Position posMM = pcP[i].position ;
4978 
4979  for (size_t j=0; j<msg->numAllAtoms; ++j ) {
4980 
4981  BigReal qmCharge = atmP[j].charge ;
4982 
4983  BigReal force = pntCharge*qmCharge*constants ;
4984 
4985  Position rVec = posMM - atmP[j].position ;
4986 
4987  force /= rVec.length2();
4988 
4989  // We accumulate the total force felt by a point charge
4990  // due to the QM charge distribution. This total force
4991  // will be applied to the point charge if it is a real one,
4992  // or will be distirbuted over MM1 and MM2 point charges, it
4993  // this is a virtual point charge.
4994  totalForce += force*rVec.unit();
4995  }
4996 
4997  if (pcP[i].type == QMPCTYPE_CLASSICAL) {
4998  // Checking pcP was not a QM atom in another region
4999  // If so the interaction is accounted in the other region
5000  if (qmAtomGroup[pcP[i].id] == 0) {
5001  // If we already ignored MM1 charges, we take all other
5002  // non-virtual charges and apply forces directly to them.
5003  resForce[pcIndx].force += totalForce;
5004  }
5005  }
5006  else {
5007  // If we are handling virtual PC, we distribute the force over
5008  // MM1 and MM2.
5009 
5010  Force mm1Force(0), mm2Force(0);
5011 
5012  // Virtual PC are bound to MM2.
5013  int mm2Indx = pcP[i].bountToIndx;
5014  // MM2 charges are bound to MM1.
5015  int mm1Indx = pcP[mm2Indx].bountToIndx;
5016 
5017  Real Cq = pcP[i].dist;
5018 
5019  mm1Force = (1-Cq)*totalForce ;
5020  mm2Force = Cq*totalForce ;
5021 
5022  resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
5023  resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
5024  }
5025 
5026  }
5027  }
5028 
5029  // Adjusts forces from PME, canceling contributions from the QM and
5030  // direct Coulomb forces calculated here.
5031  if (msg->PMEOn) {
5032 
5033  DebugM(1,"Correcting forces and energy for PME.\n");
5034 
5035  ewaldcof = msg->PMEEwaldCoefficient;
5036  BigReal TwoBySqrtPi = 1.12837916709551;
5037  pi_ewaldcof = TwoBySqrtPi * ewaldcof;
5038 
5039  for (size_t i=0; i < msg->numQMAtoms; i++) {
5040 
5041  BigReal p_i_charge = atmP[i].charge ;
5042  Position pos_i = atmP[i].position ;
5043 
5044  for (size_t j=i+1; j < msg->numQMAtoms; j++) {
5045 
5046  BigReal p_j_charge = atmP[j].charge ;
5047 
5048  Position pos_j = atmP[j].position ;
5049 
5050  BigReal r = Vector(pos_i - pos_j).length() ;
5051 
5052  BigReal tmp_a = r * ewaldcof;
5053  BigReal tmp_b = erfc(tmp_a);
5054  BigReal corr_energy = tmp_b;
5055  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
5056 
5057 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
5058  BigReal recip_energy = (1-tmp_b)/r;
5059 
5060 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
5061  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
5062 
5063  const BigReal kq_i = p_i_charge * constants;
5064 
5065  // Final force and energy correction for this pair of atoms.
5066  BigReal energy = kq_i * p_j_charge * recip_energy ;
5067 
5068  Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
5069 
5070  // The force is *subtracted* from the total force acting on
5071  // both atoms. The sign on fixForce corrects the orientation
5072  // of the subtracted force.
5073 // DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
5074 // << std::endl);
5075 // DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
5076 // << std::endl);
5077 // DebugM(4,"Force correction: " << fixForce << std::endl);
5078  resForce[i].force -= fixForce ;
5079  resForce[j].force -= -1*fixForce;
5080 
5081  // The energy is *subtracted* from the total energy calculated here.
5082  resMsg->energyCorr -= energy;
5083  }
5084  }
5085 
5086  pcIndx = msg->numQMAtoms;
5087  // We only loop over point charges from real atoms, ignoring the ones
5088  // created to handle QM-MM bonds.
5089  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
5090 
5091  BigReal p_i_charge = pcPpme[i].charge;
5092  Position pos_i = pcPpme[i].position ;
5093 
5094  const BigReal kq_i = p_i_charge * constants;
5095 
5096  Force fixForce = 0;
5097 
5098  for (size_t j=0; j<msg->numQMAtoms; ++j ) {
5099 
5100  BigReal p_j_charge = atmP[j].charge ;
5101 
5102  Position pos_j = atmP[j].position ;
5103 
5104  BigReal r = Vector(pos_i - pos_j).length() ;
5105 
5106  BigReal tmp_a = r * ewaldcof;
5107  BigReal tmp_b = erfc(tmp_a);
5108  BigReal corr_energy = tmp_b;
5109  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
5110 
5111 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
5112  BigReal recip_energy = (1-tmp_b)/r;
5113 
5114 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
5115  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
5116 
5117  // Final force and energy correction for this pair of atoms.
5118  BigReal energy = kq_i * p_j_charge * recip_energy ;
5119 
5120  fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
5121 
5122  // The energy is *subtracted* from the total energy calculated here.
5123  resMsg->energyCorr -= energy;
5124 
5125  }
5126 
5127  // The force is *subtracted* from the total force acting on
5128  // the point charge.
5129 // DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
5130 // << std::endl);
5131 // DebugM(4,"Force correction: " << fixForce << std::endl);
5132  resForce[pcIndx].force -= kq_i*fixForce ;
5133  }
5134 
5135  }
5136 
5137  DebugM(1,"Determining virial...\n");
5138 
5139  // Calculates the virial contribution form the forces on QM atoms and
5140  // point charges calculated here.
5141  for (size_t i=0; i < msg->numQMAtoms; i++) {
5142  // virial used by NAMD is -'ve of normal convention, so reverse it!
5143  // virial[i][j] in file should be sum of -1 * f_i * r_j
5144  for ( int k=0; k<3; ++k )
5145  for ( int l=0; l<3; ++l )
5146  resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
5147  }
5148 
5149  pcIndx = msg->numQMAtoms; // Index in the real PC force array.
5150  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
5151  // virial used by NAMD is -'ve of normal convention, so reverse it!
5152  // virial[i][j] in file should be sum of -1 * f_i * r_j
5153  for ( int k=0; k<3; ++k )
5154  for ( int l=0; l<3; ++l )
5155  resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
5156  }
5157 
5158  DebugM(1,"End of QM processing. Sending result message.\n");
5159 
5160  // Send message to rank zero with results.
5161  QMProxy[0].recvQMRes(resMsg);
5162 
5163  if (msg->switching && msg->PMEOn)
5164  delete [] pcPpme;
5165 
5166  delete msg;
5167  return ;
5168 }
5169 
5170 
5171 void ComputeQMMgr::pntChrgSwitching(QMGrpCalcMsg *msg, QMAtomData *pcPpme) {
5172 
5173  // We apply a switching function to the point charges so that there is a
5174  // smooth decay of the electrostatic environment seen by the QM system.
5175 
5176  BigReal cutoff2 = msg->cutoff*msg->cutoff;
5177  BigReal swdist = msg->swdist;
5178 
5179  SortedArray<pntChrgDist> sortedDists;
5180 
5181  DebugM(1,"Initiating point charge switching and processing in rank "
5182  << CkMyPe() << "\n" ) ;
5183 
5184  QMAtomData *pcP = msg->data + msg->numAllAtoms;
5185 
5186 #ifdef DEBUG_QM
5187  Real PCScaleCharge = 0;
5188  for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
5189  PCScaleCharge += pcP[i].charge;
5190  }
5191  DebugM(4, "The initial total Point-Charge charge is " << PCScaleCharge
5192  << " before scaling type: " << msg->switchType << "\n" );
5193 #endif
5194 
5195  // If PME is on, we return an unchanged vector of charges so that a
5196  // PME correction can be calculated.
5197  if (msg->PMEOn) {
5198  for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
5199  pcPpme[i] = pcP[i];
5200  }
5201  }
5202 
5203  switch (msg->switchType) {
5204 
5205  case QMPCSCALESHIFT:
5206  {
5207 
5208  // We store all point charges so that only the furthest away
5209  // are changed in PC schemes "round" or "zero"
5210  for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
5211  sortedDists.add( pntChrgDist(i, pcP[i].dist) ) ;
5212 
5213  // https://nmr.cit.nih.gov/xplor-nih/xplorMan/node119.html
5214  // Applying X-PLOR shifting formula:
5215  // F = Qi*Qj*(C/(epsilon*r))*(1 - r^2/cutoff^2)^2
5216  Real dist2 = pcP[i].dist*pcP[i].dist ;
5217  dist2 /= cutoff2 ;
5218  Real coef = 1- dist2;
5219  pcP[i].charge *= coef*coef ;
5220  }
5221 
5222  } break;
5223 
5224  case QMPCSCALESWITCH:
5225  {
5226  for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
5227 
5228  // We store the point charges which are beiond the switching
5229  // distance.
5230  if (pcP[i].dist > swdist) {
5231  sortedDists.add( pntChrgDist(i, pcP[i].dist) ) ;
5232 
5233  // https://nmr.cit.nih.gov/xplor-nih/xplorMan/node118.html
5234  // Applying the XPLOR switching formula:
5235  // (r^2 - cutoff^2)^2*(cutoff^2 + 2*r^2 - 3*swdits^2)/(cutoff^2 - swdits^2)^3
5236  Real dist2 = pcP[i].dist*pcP[i].dist ;
5237  Real swdist2 = swdist*swdist;
5238  Real coef = (dist2 - cutoff2)*(dist2 - cutoff2) ;
5239  coef *= (cutoff2 + 2*dist2 - 3*swdist2) ;
5240  coef /= (cutoff2 - swdist2)*(cutoff2 - swdist2)*(cutoff2 - swdist2);
5241  pcP[i].charge *= coef ;
5242  }
5243  }
5244  } break;
5245  }
5246 
5247 #ifdef DEBUG_QM
5248  PCScaleCharge = 0;
5249  for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
5250  PCScaleCharge += pcP[i].charge;
5251  }
5252  DebugM(4, "The final total Point-Charge charge is " << PCScaleCharge
5253  << " after scaling.\n" );
5254 #endif
5255 
5256  DebugM(4, sortedDists.size()
5257  << " point charges were selected for point charge scheme " << msg->pcScheme << "\n" );
5258 
5259  Real totalPCCharge = 0, correction = 0;
5260  switch (msg->pcScheme) {
5261 
5262  case QMPCSCHEMEROUND:
5263  {
5264  for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
5265  totalPCCharge += pcP[i].charge;
5266  }
5267  DebugM(4, "The total Point-Charge charge is " << totalPCCharge
5268  << "\n" );
5269 
5270  if ((fabsf(roundf(totalPCCharge) - totalPCCharge) <= 0.001f) ) {
5271  DebugM(4, "Charge is already a whole number!\n" );
5272  }
5273  else {
5274  correction = roundf(totalPCCharge) -totalPCCharge ;
5275  DebugM(4, "Adding to system the charge: " << correction << "\n" );
5276  }
5277  } break;
5278 
5279  case QMPCSCHEMEZERO:
5280  {
5281  for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
5282  totalPCCharge += pcP[i].charge;
5283  }
5284  DebugM(4, "The total Point-Charge charge is " << totalPCCharge << "\n");
5285 
5286  DebugM(4, "Total QM system charge is: " << msg->charge << "\n" );
5287 
5288  correction = -1*(totalPCCharge + msg->charge);
5289  if ((fabsf(correction) <= 0.001f) ) {
5290  correction = 0;
5291  DebugM(4, "Total QM + PC charge is already zero!\n" );
5292  }
5293  else
5294  DebugM(4, "Adding a charge of " << correction << " to the system\n");
5295 
5296  } break;
5297  }
5298 
5299  if (correction != 0) {
5300 
5301  int maxi = sortedDists.size(), mini = sortedDists.size()/2;
5302  Real fraction = correction/(maxi - mini);
5303 
5304  DebugM(4, "The charge fraction added to the " << (maxi - mini)
5305  << " most distant point charges is " << fraction << "\n");
5306 
5307  for (size_t i=mini; i<maxi ; i++) {
5308 
5309  pcP[sortedDists[i].index].charge += fraction ;
5310 
5311  }
5312 
5313  #ifdef DEBUG_QM
5314  totalPCCharge = 0;
5315  for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
5316  totalPCCharge += pcP[i].charge;
5317  }
5318  DebugM(4, "The total Point-Charge charge is " << totalPCCharge
5319  << "\n");
5320  #endif
5321  }
5322 
5323 }
5324 
5325 void ComputeQMMgr::lssPrepare() {
5326 
5327  lssTotRes = 0;
5328  lssResMass = 0;
5329 
5330  DebugM (4, "Preparing LSS for " << numQMGrps << " QM groups.\n" )
5331 
5332  for (int i=0; i<numQMGrps; i++) {
5333  lssTotRes += qmLSSSize[i];
5334  }
5335 
5336  lssPos = new Position[lssTotRes];
5337 
5338  grpIDResNum.resize(numQMGrps);
5339 
5340  if (simParams->qmLSSMode == QMLSSMODECOM) {
5341 
5342  lssGrpRefMass.resize(numQMGrps);
5343 
5344  for (int i=0; i<qmLSSResSize; i++)
5345  lssResMass += qmLSSMass[i];
5346 
5347  DebugM(4, "Total residue mass is " << lssResMass << "\n" )
5348  }
5349 
5350  // Get all atom IDs of solvent QM atoms, per group.
5351  int solvResBeg = 0, refAtmBeg = 0, locIter = 0, solvIndx = 0;
5352  int *lssIndxs, *refAtmsIndxs ;
5353  Mass *lssMasses, *refAtmMasses;
5354  while (locIter < numQMGrps) {
5355  lssIndxs = qmLSSIdxs + solvResBeg;
5356 
5357  DebugM (4, "Loading atom IDs for QM group " << locIter
5358  << " with " << qmLSSSize[locIter]
5359  << " solvent molecules.\n" )
5360 
5361 
5362  switch (simParams->qmLSSMode) {
5363 
5364  case QMLSSMODECOM:
5365  {
5366  lssMasses = qmLSSMass + solvResBeg;
5367 
5368  // Loads data on QM solvent residues and their atoms.
5369  for (int i=0; i<qmLSSSize[locIter]; i++) {
5370 
5371  lssPos[solvIndx] = 0;
5372 
5373  Debug( iout << "Getting atom IDs from QM solvent molecule "
5374  << solvIndx << "\n") ;
5375  for (int j=0; j<qmLSSResSize; j++) {
5376 
5377  int atmID = lssIndxs[i*qmLSSResSize + j];
5378  Mass atmMass = lssMasses[i*qmLSSResSize + j];
5379  Debug( iout << atmID << " (" << atmMass << ") ") ;
5380 
5381  grpIDResNum[locIter].insert(atmLSSData(atmID, LSSDataStr(solvIndx,atmMass)));
5382 
5383  }
5384  solvIndx++;
5385  Debug( iout << "\n" << endi );
5386  }
5387 
5388  // Loads data on the mass of QM atoms which will be used to calculate
5389  // the COM for solvent selection.
5390  refAtmsIndxs = qmLSSRefIDs + refAtmBeg;
5391  refAtmMasses = molPtr->get_qmLSSRefMass() + refAtmBeg;
5392  for (int i=0; i<qmLSSRefSize[locIter]; i++) {
5393  lssGrpRefMass[locIter].insert(
5394  refLSSData( refAtmsIndxs[i], refAtmMasses[i] )
5395  ) ;
5396  }
5397  refAtmBeg += qmLSSRefSize[locIter] ;
5398  } break ;
5399 
5400  case QMLSSMODEDIST:
5401  {
5402  // Loads data on QM solvent residues and their atoms.
5403  for (int i=0; i<qmLSSSize[locIter]; i++) {
5404 
5405  lssPos[solvIndx] = 0;
5406 
5407  Debug( iout << "Getting atom IDs from QM solvent molecule "
5408  << solvIndx << "\n") ;
5409  for (int j=0; j<qmLSSResSize; j++) {
5410 
5411  int atmID = lssIndxs[i*qmLSSResSize + j];
5412  Debug( iout << atmID << " ") ;
5413 
5414  grpIDResNum[locIter].insert(atmLSSData(atmID, LSSDataStr(solvIndx,0)));
5415 
5416  }
5417  solvIndx++;
5418  Debug( iout << "\n" << endi );
5419  }
5420 
5421  } break ;
5422 
5423  }
5424 
5425  solvResBeg += qmLSSSize[locIter]*qmLSSResSize ;
5426  locIter++;
5427  }
5428 
5429  return ;
5430 }
5431 
5432 void ComputeQMMgr::lssUpdate(int grpIter, QMAtmVec& grpQMAtmVec,
5433  QMPCVec& grpPntChrgVec) {
5434 
5435  SortedArray<lssDistSort> solvDist;
5436 
5437  Position refCOM(0) ;
5438  Mass totMass = 0;
5439 
5440  DebugM(3, "LSS UPDATE...\n")
5441 
5442  int solvResBeg = 0 ;
5443  for (int i=0; i<grpIter; i++)
5444  solvResBeg += qmLSSSize[i] ;
5445 
5446  switch (simParams->qmLSSMode ) {
5447 
5448  case QMLSSMODECOM:
5449  {
5450  DebugM(3, "Using COM for LSS in group " << grpIter << "\n")
5451 
5452  // Determined the reference center of mass for this group.
5453  for(int i=0; i<grpQMAtmVec.size(); i++) {
5454 
5455  auto it = lssGrpRefMass[grpIter].find(grpQMAtmVec[i].id);
5456  if ( it != lssGrpRefMass[grpIter].end() ) {
5457  refCOM += grpQMAtmVec[i].position*it->second ;
5458  totMass += it->second ;
5459  }
5460  }
5461 
5462  refCOM /= totMass;
5463  DebugM ( 3, "Reference COM position: " << refCOM << "\n");
5464 
5465  // Initialize the positions of all COM of quantum residues
5466  for (int solvIter=solvResBeg; solvIter<solvResBeg+qmLSSSize[grpIter]; solvIter++) {
5467  lssPos[solvIter] = 0 ;
5468  }
5469 
5470  DebugM(3, "Calculating distance of QM solvent COM from reference COM of group.\n")
5471 
5472  // Temporary handler of lssDistSort structures, while we accumulate atoms
5473  // on their respective residues and calculate distances.
5474  std::map<int,lssDistSort > resQMDist ;
5475 
5476  // Initiates COM determination for all QM solvent molecules.
5477  for(int i=0; i<grpQMAtmVec.size(); i++) {
5478  auto it = grpIDResNum[grpIter].find(grpQMAtmVec[i].id) ;
5479  if (it != grpIDResNum[grpIter].end()) {
5480  lssPos[it->second.resIndx] += grpQMAtmVec[i].position*it->second.mass ;
5481 
5482  // tries to find a residue number
5483  auto itRes = resQMDist.find(it->second.resIndx) ;
5484  if (itRes == resQMDist.end() ) {
5485  resQMDist.insert(std::pair<int,lssDistSort >(
5486  it->second.resIndx, lssDistSort(QMLSSQMRES, -1) ));
5487  }
5488 
5489  // For each classical residue ID, we compile a list of atom IDs which
5490  // are found among point charges
5491 // resQMDist[it->second.resIndx].idVec.push_back(grpQMAtmVec[i].id) ;
5492 // resQMDist[it->second.resIndx].indxVec.push_back(i) ;
5493 
5494  resQMDist[it->second.resIndx].idIndx.add(idIndxStr(grpQMAtmVec[i].id,i)) ;
5495  }
5496  }
5497 
5498  DebugM(3, "QM Solvent molecules " << solvResBeg
5499  << " to " << solvResBeg+qmLSSSize[grpIter] << "\n")
5500 
5501 
5502  for (int solvIter=solvResBeg; solvIter<solvResBeg+qmLSSSize[grpIter]; solvIter++) {
5503  Real dist = Vector((lssPos[solvIter]/lssResMass) - refCOM).length() ;
5504  resQMDist[solvIter].dist = dist ;
5505  solvDist.add(resQMDist[solvIter]);
5506  }
5507 
5508  #ifdef DEBUG_QM
5509  DebugM(3, "We loaded the following QM solvent residues and distances:\n")
5510  for (int i=0; i<solvDist.size(); i++) {
5511  iout << i << ") type: " << solvDist[i].type
5512  << " dist " << solvDist[i].dist
5513  << " IDs: " ;
5514  for (int j=0; j<solvDist[i].idIndx.size(); j++)
5515  iout << solvDist[i].idIndx[j].ID << " " ;
5516  iout << "\n" << endi;
5517  }
5518  #endif
5519 
5520  // Get Center Of Mass distance for all (whole) solvent molecules among
5521  // *point charges*, comparing against the non-solvent QM atoms.
5522 
5523  // This will count how many PCs are associated with each solvent residue,
5524  // and associate their atom ID with that residue.
5525  std::map<int,lssDistSort > resPCSize ;
5526 
5527  for(int i=0; i<grpPntChrgVec.size(); i++) {
5528 
5529  // This maps aomt IDs with a global classical residue ID.
5530  auto it = molPtr->get_qmMMSolv().find(grpPntChrgVec[i].id) ;
5531  if (it != molPtr->get_qmMMSolv().end()) {
5532 
5533  // tries to find a residue number
5534  auto itRes = resPCSize.find(it->second) ;
5535  if (itRes == resPCSize.end() ) {
5536  resPCSize.insert(std::pair<int,lssDistSort >(
5537  it->second, lssDistSort(QMLSSCLASSICALRES, -1) ));
5538  }
5539 
5540  // For each classical residue ID, we compile a list of atom IDs which
5541  // are found among point charges
5542 // resPCSize[it->second].idVec.push_back(grpPntChrgVec[i].id) ;
5543 // resPCSize[it->second].indxVec.push_back(i) ;
5544 
5545  resPCSize[it->second].idIndx.add(idIndxStr(grpPntChrgVec[i].id,i)) ;
5546  }
5547  }
5548 
5549  // Now we check which classical solvent molecules are complete,
5550  // and compute the COM for the complete ones, while ignoring the
5551  // incomplete ones.
5552  for (auto it=resPCSize.begin(); it!=resPCSize.end(); it++) {
5553 
5554  if (it->second.idIndx.size() == qmLSSResSize) {
5555 
5556  Position currCOM(0);
5557  Mass totalMass = 0;
5558 
5559  // iout << "Found complete classical residue " << it->first << "\n" << endi;
5560 
5561  for (int i=0; i<it->second.idIndx.size(); i++) {
5562  // iout << "AtomID " << it->second.idIndx[i] .ID
5563  // << " indxVec " << it->second.idIndx[i].indx
5564  // << " mass " << grpPntChrgVec[it->second.idIndx[i].indx].mass
5565  // << " position " << grpPntChrgVec[it->second.idIndx[i].indx].position << "\n" << endi;
5566  currCOM += grpPntChrgVec[it->second.idIndx[i].indx].position*grpPntChrgVec[it->second.idIndx[i].indx].mass;
5567  totalMass += grpPntChrgVec[it->second.idIndx[i].indx].mass;
5568  }
5569  currCOM /= totalMass;
5570 
5571  Real currDist = Vector(currCOM - refCOM).length() ;
5572 
5573  it->second.dist = currDist ;
5574 
5575  solvDist.add(it->second) ;
5576  }
5577 
5578  }
5579  } break ;
5580 
5581  case QMLSSMODEDIST:
5582  {
5583  DebugM(3, "Using minimal distances for LSS in group " << grpIter << "\n")
5584 
5585  DebugM(3, "QM Solvent molecules " << solvResBeg
5586  << " to " << solvResBeg+qmLSSSize[grpIter] << "\n")
5587 
5588  // List of QM indices which will be used as reference for distance calculation.
5589  ResizeArray<int> qmRefIndx ;
5590 
5591  // Temporary handler of lssDistSort structures, while we accumulate atoms
5592  // on their respective residues and calculate distances.
5593  std::map<int,lssDistSort > resQMDist ;
5594 
5595  // Initiates COM determination for all QM solvent molecules.
5596  for(int i=0; i<grpQMAtmVec.size(); i++) {
5597 
5598  auto it = grpIDResNum[grpIter].find(grpQMAtmVec[i].id) ;
5599  if (it != grpIDResNum[grpIter].end()) {
5600 
5601  // tries to find a residue number
5602  auto itRes = resQMDist.find(it->second.resIndx) ;
5603  if (itRes == resQMDist.end() ) {
5604  resQMDist.insert(std::pair<int,lssDistSort >(
5605  it->second.resIndx, lssDistSort(QMLSSQMRES, -1) ));
5606  }
5607 
5608  // For each classical residue ID, we compile a list of atom IDs which
5609  // are found among point charges
5610 // resQMDist[it->second.resIndx].idVec.push_back(grpQMAtmVec[i].id) ;
5611 // resQMDist[it->second.resIndx].indxVec.push_back(i) ;
5612 
5613  resQMDist[it->second.resIndx].idIndx.add(idIndxStr(grpQMAtmVec[i].id,i)) ;
5614 
5615  }
5616  else {
5617  qmRefIndx.add(i) ;
5618  }
5619  }
5620 
5621  // We now calculate the shortest distance between a QM solvent residue
5622  // and any non-solvent QM atom.
5623  for (auto it=resQMDist.begin(); it != resQMDist.end(); it++) {
5624 
5625  // We prime the residue distance with the first non-solvent
5626  // QM atom.
5627  it->second.dist = Vector(
5628  grpQMAtmVec[it->second.idIndx[0].indx].position -
5629  grpQMAtmVec[qmRefIndx[0]].position
5630  ).length() ;
5631 
5632  for (int i=0; i<it->second.idIndx.size(); i++) {
5633 
5634  for(int j=0; j<qmRefIndx.size(); j++) {
5635  Real currDist = Vector(
5636  grpQMAtmVec[it->second.idIndx[i].indx].position -
5637  grpQMAtmVec[qmRefIndx[j]].position
5638  ).length() ;
5639 
5640  if (currDist < it->second.dist)
5641  it->second.dist = currDist;
5642  }
5643  }
5644 
5645  // After determining the distance of this QM solvent residue,
5646  // we add it to the sorted list.
5647  solvDist.add(it->second) ;
5648  }
5649 
5650 
5651  #ifdef DEBUG_QM
5652  DebugM(3, "We loaded the following QM solvent residues and distances:\n")
5653  for (int i=0; i<solvDist.size(); i++) {
5654  iout << i << ") type: " << solvDist[i].type
5655  << " dist " << solvDist[i].dist
5656  << " IDs: " ;
5657  for (int j=0; j<solvDist[i].idIndx.size(); j++)
5658  iout << solvDist[i].idIndx[j].ID << " " ;
5659  iout << "\n" << endi;
5660  }
5661  #endif
5662 
5663  // Get shortest distance for all (whole) solvent molecules among
5664  // *point charges*, comparing against the non-solvent QM atoms.
5665 
5666  // This will count how many PCs are associated with each solvent residue,
5667  // and associate their atom ID with that residue.
5668  std::map<int,lssDistSort > resPCSize ;
5669 
5670  for(int i=0; i<grpPntChrgVec.size(); i++) {
5671 
5672  // This maps aomt IDs with a global classical residue ID.
5673  auto it = molPtr->get_qmMMSolv().find(grpPntChrgVec[i].id) ;
5674  if (it != molPtr->get_qmMMSolv().end()) {
5675 
5676  // tries to find a residue number
5677  auto itRes = resPCSize.find(it->second) ;
5678  if (itRes == resPCSize.end() ) {
5679  resPCSize.insert(std::pair<int,lssDistSort >(
5680  it->second, lssDistSort(QMLSSCLASSICALRES, -1) ));
5681  }
5682 
5683  // For each classical residue ID, we compile a list of atom IDs which
5684  // are found among point charges
5685 // resPCSize[it->second].idVec.push_back(grpPntChrgVec[i].id) ;
5686 // resPCSize[it->second].indxVec.push_back(i) ;
5687 
5688  resPCSize[it->second].idIndx.add(idIndxStr(grpPntChrgVec[i].id,i)) ;
5689  }
5690  }
5691 
5692  // Now we check which classical solvent molecules are complete,
5693  // and compute the COM for the complete ones, while ignoring the
5694  // incomplete ones.
5695  for (auto it=resPCSize.begin(); it!=resPCSize.end(); it++) {
5696 
5697  if (it->second.idIndx.size() == qmLSSResSize) {
5698 
5699  // We prime the residue distance with the first non-solvent
5700  // QM atom.
5701  it->second.dist = Vector(
5702  grpPntChrgVec[it->second.idIndx[0].indx].position -
5703  grpQMAtmVec[qmRefIndx[0]].position
5704  ).length() ;
5705 
5706  for (int i=0; i<it->second.idIndx.size(); i++) {
5707 
5708  for(int j=0; j<qmRefIndx.size(); j++) {
5709  Real currDist = Vector(
5710  grpPntChrgVec[it->second.idIndx[i].indx].position -
5711  grpQMAtmVec[qmRefIndx[j]].position
5712  ).length() ;
5713 
5714  if (currDist < it->second.dist)
5715  it->second.dist = currDist;
5716  }
5717  }
5718 
5719  solvDist.add(it->second) ;
5720  }
5721 
5722  }
5723 
5724  } break ;
5725 
5726  } // End switch
5727 
5728  #ifdef DEBUG_QM
5729  DebugM(3, "Final selection of solvent residues and distances:\n")
5730  for (int i=0; i<qmLSSSize[grpIter]; i++) {
5731  std::string typeS ;
5732  if (solvDist[i].type != QMLSSCLASSICALRES )
5733  typeS = "QM" ;
5734  else
5735  typeS = "Classical";
5736  iout << i << ") type: " << typeS
5737  << " dist " << solvDist[i].dist
5738  << " IDs: " ;
5739  for (int j=0; j<solvDist[i].idIndx.size(); j++)
5740  iout << solvDist[i].idIndx[j].ID << " " ;
5741  iout << "\n" << endi;
5742  }
5743  #endif
5744 
5745  // Compare COM distances of QM and Classical solvent molecules, creating a
5746  // substitution list, atmID by atmID.
5747 
5748  DebugM(3, "Determining residues to be swaped...\n")
5749 
5750  ResizeArray<lssDistSort> nearMM, farQM ;
5751 
5752  for (int resIter = 0; resIter < qmLSSSize[grpIter] ; resIter++) {
5753  if (solvDist[resIter].type == QMLSSCLASSICALRES) {
5754  nearMM.add(solvDist[resIter]);
5755  }
5756  }
5757 
5758  for (int resIter=qmLSSSize[grpIter]; resIter<solvDist.size(); resIter++) {
5759  if (solvDist[resIter].type == QMLSSQMRES) {
5760  farQM.add(solvDist[resIter]);
5761  }
5762 
5763  if (farQM.size() == nearMM.size()) break;
5764  }
5765 
5766  if (farQM.size() != nearMM.size())
5767  NAMD_die("Could not find complementing residues to be swapped in LSS.") ;
5768 
5769  #ifdef DEBUG_QM
5770  DebugM(3, "Removing the following QM residues:\n")
5771  for (int i=0; i<farQM.size();i++) {
5772  std::string typeS ;
5773  if (farQM[i].type != QMLSSCLASSICALRES )
5774  typeS = "QM" ;
5775  else
5776  typeS = "Classical";
5777  iout << i << ") type: " << typeS
5778  << " dist " << farQM[i].dist
5779  << " IDs: " ;
5780  for (int j=0; j<farQM[i].idIndx.size(); j++)
5781  iout << farQM[i].idIndx[j].ID << " " ;
5782  iout << "\n" << endi;
5783  }
5784 
5785  DebugM(3, "Replacing with the following Classical residues:\n")
5786  for (int i=0; i<nearMM.size();i++) {
5787  std::string typeS ;
5788  if (nearMM[i].type != QMLSSCLASSICALRES )
5789  typeS = "QM" ;
5790  else
5791  typeS = "Classical";
5792  iout << i << ") type: " << typeS
5793  << " dist " << nearMM[i].dist
5794  << " IDs: " ;
5795  for (int j=0; j<nearMM[i].idIndx.size(); j++)
5796  iout << nearMM[i].idIndx[j].ID << " " ;
5797  iout << "\n" << endi;
5798  }
5799 
5800  DebugM(3, "Building substitution array...\n")
5801  #endif
5802 
5803  iout << iINFO << "LSS is swapping " << farQM.size() << " solvent residues in QM group "
5804  << grpIter << "\n" << endi;
5805 
5806  // Now we build the array which will be sent to all nodes with force results from
5807  // this step.
5808  // Atom reassignment will be done in the next step, and will use this data.
5809  for (int i=0; i<farQM.size();i++) {
5810 
5811  for(int j=0; j<qmLSSResSize; j++) {
5812 
5813  int qIndx= farQM[i].idIndx[j].indx;
5814  int mIndx= nearMM[i].idIndx[j].indx;
5815 
5816  subsArray.add( LSSSubsDat( grpQMAtmVec[qIndx].id,
5817  grpPntChrgVec[mIndx].id,
5818  grpPntChrgVec[mIndx].vdwType,
5819  grpPntChrgVec[mIndx].charge ) );
5820 
5821  subsArray.add( LSSSubsDat( grpPntChrgVec[mIndx].id,
5822  grpQMAtmVec[qIndx].id,
5823  grpQMAtmVec[qIndx].vdwType,
5824  grpQMAtmVec[qIndx].charge ) );
5825  }
5826 
5827  }
5828 
5829  #ifdef DEBUG_QM
5830  for(int i=0; i<subsArray.size() ;i++) {
5831  DebugM(3, CkMyPe() << ") Storing LSS atom " << subsArray[i].origID
5832  << " Which will become " << subsArray[i].newID
5833  << " - " << subsArray[i].newVdWType
5834  << " - " << subsArray[i].newCharge << "\n" );
5835  }
5836  #endif
5837 
5838  return;
5839 }
5840 
5841 void ComputeQMMgr::calcCSMD(int grpIndx, int numQMAtoms, const QMAtomData *atmP, Force *cSMDForces) {
5842 
5843  // For each Conditional SMD instance, we update the guide particle
5844  // and calculate and apply the respective force.
5845  for ( int i=0; i < cSMDindxLen[grpIndx]; i++ ) {
5846 
5847  // Defite target distance
5848  Real tdist;
5849  Force cSMDforce(0);
5850 
5851  // cSMD instance index
5852  int cSMDit = cSMDindex[grpIndx][i] ;
5853  AtomID src = -1, trgt = -1;
5854 
5855  // Finds local indices for the atoms involved in this cSMD instance
5856  for (size_t indx=0; indx < numQMAtoms; ++indx) {
5857 
5858  if ( cSMDpairs[cSMDit][0] == atmP[indx].id )
5859  src = indx ;
5860 
5861  if ( cSMDpairs[cSMDit][1] == atmP[indx].id )
5862  trgt = indx ;
5863 
5864  if (src != -1 && trgt != -1)
5865  break;
5866  }
5867 
5868  // Finds the unit vector in the direction of the target
5869  Vector trgtDir = atmP[trgt].position - atmP[src].position ;
5870  tdist = trgtDir.length();
5871  trgtDir /= tdist;
5872 
5873  // If this is the first time step, set up the guide particle.
5874  if (cSMDInitGuides < cSMDnumInstances) {
5875  cSMDguides[cSMDit] = atmP[src].position;
5876  cSMDInitGuides++;
5877  }
5878 
5879  // If the target is further away than "cutoff", advance the
5880  // guide particle and apply force.
5881  if (tdist > cSMDcoffs[cSMDit]) {
5882 
5883  cSMDguides[cSMDit] += trgtDir*cSMDVels[cSMDit];
5884 
5885  Vector guideDir = cSMDguides[cSMDit] - atmP[src].position ;
5886 
5887  // Calculate force in the direction of the guide particle.
5888  guideDir *= cSMDKs[cSMDit] ;
5889  cSMDforce = guideDir ;
5890 
5891  }
5892  else {
5893  // If we pulled the Src towards Trgt closer than Cutoff,
5894  // we reset the guide particle to the position of Src.
5895  cSMDguides[cSMDit] = atmP[src].position;
5896  }
5897 
5898  DebugM(4, cSMDit << ") Center at " << cSMDguides[cSMDit] << " for target distance " <<
5899  tdist << " and direction " << trgtDir <<
5900  "." << std::endl);
5901 
5902  // We now save the force in a vector of length "cSMDnumInstances".
5903  cSMDForces[cSMDit] = cSMDforce;
5904 
5905  iout << iINFO << "Calculated cSMD force vector " << cSMDforce
5906  << " (atom " << cSMDpairs[cSMDit][0] << " to " << cSMDpairs[cSMDit][1] << ")\n" << endi;
5907  }
5908 }
5909 
5910 
5911 
5912 #ifdef DEBUG_QM
5913 
5914 void ComputeQMMgr::Write_PDB(std::string Filename, const QMGrpCalcMsg *dataMsg)
5915 {
5916  std::ofstream OutputTmpPDB ;
5917 
5918  try
5919  {
5920 
5921  OutputTmpPDB.open( Filename.c_str(), std::fstream::out );
5922 
5923  OutputTmpPDB << "REMARK Information used by NAMD to create the QM simulation." << std::endl ;
5924  OutputTmpPDB << "REMARK Occupancy: bount to index; Beta: System ID; " << std::endl ;
5925 
5926  const QMAtomData *dataP = dataMsg->data;
5927 
5928  for ( int i=0 ; i < dataMsg->numAllAtoms + dataMsg->numAllPntChrgs ; i++ )
5929  {
5930  // PDB format: http://www.wwpdb.org/documentation/format33/sect9.html
5931 
5932  // Exception: The field resName was changed from the coluns 18-20 to the columns 18-21
5933  // This allows the usage of protonated amino acids and other molecules in the PDB file.
5934 
5935  // Exception2: The field extraInfo was changed from the coluns 67 - 76 to the columns 67 - 72
5936  // This allows the usage of segments in PDB files, necessary for the propper usage of NAMD.
5937 
5938  std::string name(" x ");
5939  std::string resName (" uk");
5940  std::string chainName("X");
5941  std::string element("") ;
5942  if (i < dataMsg->numAllAtoms ) {
5943  name = dataP[i].element;
5944  chainName = "q" ;
5945  element = dataP[i].element;
5946  if (dataP[i].type == QMATOMTYPE_QM)
5947  resName = " qm " ;
5948  else if (dataP[i].type == QMATOMTYPE_DUMMY)
5949  resName = " dm " ;
5950  }
5951  else {
5952  chainName = "c" ;
5953  if (dataP[i].type == QMPCTYPE_CLASSICAL)
5954  resName = " pc ";
5955  else if (dataP[i].type == QMPCTYPE_IGNORE)
5956  resName = "ipc ";
5957  else if (dataP[i].type == QMPCTYPE_EXTRA)
5958  resName = "vpc ";
5959  }
5960 
5961  OutputTmpPDB << "ATOM " ; // ATOM 1 - 6
5962 
5963  OutputTmpPDB.width(5) ; // serial 7 - 11
5964  OutputTmpPDB.right ;
5965  OutputTmpPDB << i ;
5966 
5967  OutputTmpPDB << " " ; // Spacing
5968 
5969  OutputTmpPDB.width(4) ; // name 13 - 16
5970  OutputTmpPDB << name ;
5971 
5972  OutputTmpPDB << " " ; // altLoc 17
5973 
5974  OutputTmpPDB.width(4) ; // resName 18 - 21
5975  OutputTmpPDB << resName ;
5976 
5977  OutputTmpPDB.width(1) ; // chainID 22
5978  OutputTmpPDB << chainName ;
5979 
5980  OutputTmpPDB.width(4) ; // Residue Index 23 - 26
5981  OutputTmpPDB << i ;
5982 
5983  OutputTmpPDB << " " ; // iCode 27
5984 
5985  OutputTmpPDB << " " ; // Spacing
5986 
5987  OutputTmpPDB.width(8) ; // x 31 - 38
5988  OutputTmpPDB.right ;
5989  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
5990  OutputTmpPDB.precision(3) ;
5991  OutputTmpPDB << dataP[i].position.x ;
5992 
5993  OutputTmpPDB.width(8) ; // y 39 - 46
5994  OutputTmpPDB.right ;
5995  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
5996  OutputTmpPDB.precision(3) ;
5997  OutputTmpPDB << dataP[i].position.y ;
5998 
5999  OutputTmpPDB.width(8) ; // z 47 - 54
6000  OutputTmpPDB.right ;
6001  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
6002  OutputTmpPDB.precision(3) ;
6003  OutputTmpPDB << dataP[i].position.z ;
6004 
6005  OutputTmpPDB.width(6) ; // occupancy 55 - 60
6006  OutputTmpPDB.right ;
6007  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
6008  OutputTmpPDB.precision(2) ;
6009  OutputTmpPDB << dataP[i].bountToIndx ;
6010 
6011  OutputTmpPDB.width(6) ; // tempFactor/Beta 61 - 66
6012  OutputTmpPDB.right ;
6013  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
6014  OutputTmpPDB.precision(2) ;
6015  OutputTmpPDB << dataP[i].id ;
6016 
6017  OutputTmpPDB.width(6) ; // extra information not originaly on PDB format 67 - 72
6018  OutputTmpPDB << " " ;
6019 
6020  OutputTmpPDB.width(4) ; // segment information from NAMD, not originaly on PDB format 72 - 76
6021  OutputTmpPDB.left ;
6022  OutputTmpPDB << "QM ";
6023 
6024  OutputTmpPDB.width(2) ; // element 77 - 78
6025  OutputTmpPDB.right ;
6026  OutputTmpPDB << element ;
6027 
6028  OutputTmpPDB.width(2) ; // charge 77 - 78
6029  OutputTmpPDB.right ;
6030  OutputTmpPDB << dataP[i].charge ;
6031 
6032  OutputTmpPDB << std::endl;
6033 
6034  }
6035 
6036  OutputTmpPDB << "END" << std::endl;
6037 
6038  OutputTmpPDB.close();
6039  }
6040  catch (...)
6041  {
6042  iout << iERROR << "Generic exception at QM write PBD." << endi ;
6043  NAMD_die("PDB write error");
6044  throw "Generic exception!" ;
6045  }
6046  return ;
6047 }
6048 
6049 void ComputeQMMgr::Write_PDB(std::string Filename, const QMCoordMsg *dataMsg)
6050 {
6051  std::ofstream OutputTmpPDB ;
6052 
6053  try
6054  {
6055 
6056  OutputTmpPDB.open( Filename.c_str(), std::fstream::out );
6057 
6058  OutputTmpPDB << "REMARK Information used by NAMD to create the QM simulation." << std::endl ;
6059  OutputTmpPDB << "REMARK Occupancy: bount to index; Beta: System ID; " << std::endl ;
6060 
6061  const ComputeQMAtom *dataP = dataMsg->coord;
6062 
6063  for ( int i=0 ; i < dataMsg->numAtoms; i++ )
6064  {
6065  // PDB format: http://www.wwpdb.org/documentation/format33/sect9.html
6066 
6067  // Exception: The field resName was changed from the coluns 18-20 to the columns 18-21
6068  // This allows the usage of protonated amino acids and other molecules in the PDB file.
6069 
6070  // Exception2: The field extraInfo was changed from the coluns 67 - 76 to the columns 67 - 72
6071  // This allows the usage of segments in PDB files, necessary for the propper usage of NAMD.
6072 
6073  std::string name(" x ");
6074  std::string resName (" atm");
6075  std::string chainName("X");
6076  std::string element("") ;
6077 
6078  OutputTmpPDB << "ATOM " ; // ATOM 1 - 6
6079 
6080  OutputTmpPDB.width(5) ; // serial 7 - 11
6081  OutputTmpPDB.right ;
6082  OutputTmpPDB << i ;
6083 
6084  OutputTmpPDB << " " ; // Spacing
6085 
6086  OutputTmpPDB.width(4) ; // name 13 - 16
6087  OutputTmpPDB << name ;
6088 
6089  OutputTmpPDB << " " ; // altLoc 17
6090 
6091  OutputTmpPDB.width(4) ; // resName 18 - 21
6092  OutputTmpPDB << resName ;
6093 
6094  OutputTmpPDB.width(1) ; // chainID 22
6095  OutputTmpPDB << chainName ;
6096 
6097  OutputTmpPDB.width(4) ; // Residue Index 23 - 26
6098  OutputTmpPDB << i ;
6099 
6100  OutputTmpPDB << " " ; // iCode 27
6101 
6102  OutputTmpPDB << " " ; // Spacing
6103 
6104  OutputTmpPDB.width(8) ; // x 31 - 38
6105  OutputTmpPDB.right ;
6106  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
6107  OutputTmpPDB.precision(3) ;
6108  OutputTmpPDB << dataP[i].position.x ;
6109 
6110  OutputTmpPDB.width(8) ; // y 39 - 46
6111  OutputTmpPDB.right ;
6112  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
6113  OutputTmpPDB.precision(3) ;
6114  OutputTmpPDB << dataP[i].position.y ;
6115 
6116  OutputTmpPDB.width(8) ; // z 47 - 54
6117  OutputTmpPDB.right ;
6118  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
6119  OutputTmpPDB.precision(3) ;
6120  OutputTmpPDB << dataP[i].position.z ;
6121 
6122  OutputTmpPDB.width(6) ; // occupancy 55 - 60
6123  OutputTmpPDB.right ;
6124  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
6125  OutputTmpPDB.precision(2) ;
6126  OutputTmpPDB << dataP[i].qmGrpID ;
6127 
6128  OutputTmpPDB.width(6) ; // tempFactor/Beta 61 - 66
6129  OutputTmpPDB.right ;
6130  OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
6131  OutputTmpPDB.precision(2) ;
6132  OutputTmpPDB << dataP[i].id ;
6133 
6134  OutputTmpPDB.width(6) ; // extra information not originaly on PDB format 67 - 72
6135  OutputTmpPDB << " " ;
6136 
6137  OutputTmpPDB.width(4) ; // segment information from NAMD, not originaly on PDB format 72 - 76
6138  OutputTmpPDB.left ;
6139  OutputTmpPDB << "QM ";
6140 
6141  OutputTmpPDB.width(2) ; // element 77 - 78
6142  OutputTmpPDB.right ;
6143  OutputTmpPDB << element ;
6144 
6145  OutputTmpPDB.width(2) ; // charge 77 - 78
6146  OutputTmpPDB.right ;
6147  OutputTmpPDB << dataP[i].charge ;
6148 
6149  OutputTmpPDB << std::endl;
6150 
6151  }
6152 
6153  OutputTmpPDB << "END" << std::endl;
6154 
6155  OutputTmpPDB.close();
6156  }
6157  catch (...)
6158  {
6159  iout << iERROR << "Generic exception at QM write PBD." << endi ;
6160  NAMD_die("PDB write error");
6161  throw "Generic exception!" ;
6162  }
6163  return ;
6164 }
6165 
6166 #endif
6167 
6168 #include "ComputeQMMgr.def.h"
6169 
static Node * Object()
Definition: Node.h:86
#define QMPCSCALESHIFT
Definition: ComputeQM.C:62
lssDistSort(const lssDistSort &ref)
Definition: ComputeQM.C:256
int bountToIndx
Definition: ComputeQM.C:299
int replace
Definition: ComputeQM.h:102
Bool get_noPC()
Definition: Molecule.h:857
int * get_qmLSSRefSize()
Definition: Molecule.h:851
BigReal PMEEwaldCoefficient
Definition: ComputeQM.C:336
#define Debug(x)
Definition: Debug.h:58
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:81
QMForce * force
Definition: ComputeQM.C:181
std::map< int, Mass > LSSRefMap
Definition: ComputeQM.C:396
ComputeQMPntChrg(Position posInit, float chrgInit, int idInit, Real qmInit, int hiInit, Real newDist, Mass newM, int newvdwType)
Definition: ComputeQM.C:136
int timestep
Definition: ComputeQM.C:107
#define QMLSSMODECOM
Definition: Molecule.h:127
const int * get_qmMMNumTargs()
Definition: Molecule.h:843
char baseDir[256]
Definition: ComputeQM.C:338
#define QMRATIOTYPE
Definition: ComputeQM.C:56
char execPath[256]
Definition: ComputeQM.C:338
#define QMSCHEMEZ2
Definition: Molecule.h:140
void NAMD_err(const char *err_msg)
Definition: common.C:106
Real * get_qmGrpMult()
Definition: Molecule.h:832
QMAtomData * data
Definition: ComputeQM.C:339
void recvPartQM(QMCoordMsg *)
Definition: ComputeQM.C:818
int get_qmNumGrps() const
Definition: Molecule.h:828
void setCompute(ComputeQM *c)
Definition: ComputeQM.C:406
int numLssDat
Definition: ComputeQM.C:190
Mass mass
Definition: ComputeQM.C:360
ComputeQMAtom(Position posInit, float chrgInit, int idInit, Real qmInit, int hiInit, int newvdwType)
Definition: ComputeQM.C:74
Real multiplicity
Definition: ComputeQM.C:327
SortedArray< LSSSubsDat > & lssSubs(ComputeQMMgr *mgr)
Definition: ComputeQM.C:595
int ComputeID
Definition: NamdTypes.h:183
idIndxStr & operator=(const idIndxStr &ref)
Definition: ComputeQM.C:221
BigReal energy
Definition: ComputeQM.C:187
int AtomID
Definition: NamdTypes.h:29
const int * get_qmGrpNumBonds()
Definition: Molecule.h:836
const int * get_qmAtmIndx()
Definition: Molecule.h:824
Vector Force
Definition: NamdTypes.h:26
#define QMSCHEMERCD
Definition: Molecule.h:138
BigReal constants
Definition: ComputeQM.C:328
#define PCMODEUPDATEPOS
Definition: ComputeQM.C:100
std::pair< Position, BigReal > PosDistPair
Definition: ComputeQM.C:1290
static PatchMap * Object()
Definition: PatchMap.h:27
#define QMPCSCHEMEROUND
Definition: ComputeQM.C:59
virtual ~ComputeQM()
Definition: ComputeQM.C:609
bool operator!=(const idIndxStr &ref) const
Definition: ComputeQM.C:230
int get_cSMDnumInst()
Definition: Molecule.h:869
int find(const Elem &e) const
Definition: ResizeArray.h:137
int add(const Elem &elem)
Definition: SortedArray.h:55
float charge
Definition: ComputeQM.C:297
int sourceNode
Definition: ComputeQM.C:169
Definition: Vector.h:64
bool prepProcOn
Definition: ComputeQM.C:330
int insert(const Elem &elem, int index)
Definition: ResizeArray.h:109
SimParameters * simParameters
Definition: Node.h:178
BigReal nonbondedScaling
ComputeHomePatchList patchList
int get_numQMAtoms()
Definition: Molecule.h:826
int numForces
Definition: ComputeQM.C:189
Mass * get_qmLSSMass()
Definition: Molecule.h:847
int qmInd
Definition: ComputeQM.C:348
float charge
Definition: ComputeQM.C:67
#define QMATOMTYPE_DUMMY
Definition: ComputeQM.C:195
int * pcIndxs
Definition: ComputeQM.C:111
float Real
Definition: common.h:109
#define COULOMB
Definition: common.h:46
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
int numAllAtoms
Definition: ComputeQM.C:324
const int *const * get_qmMMBondedIndx()
Definition: Molecule.h:839
#define QMPCTYPE_CLASSICAL
Definition: ComputeQM.C:198
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
BigReal z
Definition: Vector.h:66
#define X
Definition: msm_defn.h:29
int replace
Definition: NamdTypes.h:201
char qmPrepProc[NAMD_FILENAME_BUFFER_SIZE]
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:82
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
Real dist
Definition: ComputeQM.C:302
bool operator<(const idIndxStr &ref)
Definition: ComputeQM.C:238
bool operator==(const idIndxStr &ref) const
Definition: ComputeQM.C:227
idIndxStr(const idIndxStr &ref)
Definition: ComputeQM.C:216
Position position
Definition: ComputeQM.C:296
dummyData(Position newPos, int newQmInd, int newIndx)
Definition: ComputeQM.C:351
#define QMFormatMOPAC
Definition: Molecule.h:130
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
#define iout
Definition: InfoStream.h:51
int const * get_cSMDindxLen()
Definition: Molecule.h:870
SortedArray< idIndxStr > idIndx
Definition: ComputeQM.C:246
void storeQMRes(QMGrpResMsg *)
Definition: ComputeQM.C:2422
void recvPntChrg(QMPntChrgMsg *)
Definition: ComputeQM.C:1910
Bool get_qmcSMD()
Definition: Molecule.h:868
BigReal length(void) const
Definition: Vector.h:169
int * get_qmLSSRefIDs()
Definition: Molecule.h:850
SortedArray< LSSSubsDat > & get_subsArray()
Definition: ComputeQM.C:428
int * get_qmLSSIdxs()
Definition: Molecule.h:846
std::map< int, int > & get_qmMMSolv()
Definition: Molecule.h:853
const char *const * get_qmDummyElement()
Definition: Molecule.h:840
Position position
Definition: ComputeQM.C:126
void doWork()
Definition: ComputeQM.C:675
char qmExecPath[NAMD_FILENAME_BUFFER_SIZE]
int switchType
Definition: ComputeQM.C:333
#define QMCHRGCHELPG
Definition: Molecule.h:135
char outputFilename[NAMD_FILENAME_BUFFER_SIZE]
BigReal switchingDist
int open_dcd_write(const char *dcdname)
Definition: dcdlib.C:662
#define QMFormatORCA
Definition: Molecule.h:129
int get_qmTotCustPCs()
Definition: Molecule.h:864
Real * get_qmGrpChrg()
Definition: Molecule.h:831
ComputeQMAtom * coord
Definition: ComputeQM.C:110
#define QMFormatUSR
Definition: Molecule.h:131
const int *const * get_qmGrpBonds()
Definition: Molecule.h:838
int const *const * get_cSMDpairs()
Definition: Molecule.h:872
int insert(const Elem &elem)
#define DCD_FILEEXISTS
Definition: dcdlib.h:52
char qmSecProc[NAMD_FILENAME_BUFFER_SIZE]
bool switching
Definition: ComputeQM.C:332
Charge charge
Definition: NamdTypes.h:54
#define QMSCHEMEZ3
Definition: Molecule.h:141
lssDistSort(int newType, Real newDist)
Definition: ComputeQM.C:251
#define QMSCHEMEZ1
Definition: Molecule.h:139
static char * QMETITLE(int X)
Definition: ComputeQM.C:384
Real * get_qmMeQMGrp()
Definition: Molecule.h:860
ResizeArrayIter< T > end(void) const
pntChrgDist(int newIndx, Real newDist)
Definition: ComputeQM.C:118
const Real * get_qmAtomGroup() const
Definition: Molecule.h:820
#define QMSCHEMECS
Definition: Molecule.h:137
int get_qmMeNumBonds()
Definition: Molecule.h:858
int indx
Definition: ComputeQM.C:206
int numForces
Definition: ComputeQM.C:180
Real * get_qmGrpID()
Definition: Molecule.h:830
BigReal energyOrig
Definition: ComputeQM.C:177
lssDistSort & operator=(const lssDistSort &ref)
Definition: ComputeQM.C:264
std::pair< int, Mass > refLSSData
Definition: ComputeQM.C:395
void procQMRes()
Definition: ComputeQM.C:2485
int operator<(const pntChrgDist &v)
Definition: ComputeQM.C:122
int numQMAtoms
Definition: ComputeQM.C:323
void sort(void)
Definition: SortedArray.h:66
int id
Definition: ComputeQM.h:106
int qmAtmChrgMode
Definition: ComputeQM.C:337
bool operator<(const lssDistSort &ref)
Definition: ComputeQM.C:291
ResizeArray< ComputePme * > & getComputes(ComputePmeMgr *mgr)
Definition: ComputePme.C:593
#define PCMODEUPDATESEL
Definition: ComputeQM.C:99
int * get_qmMeMMindx()
Definition: Molecule.h:859
Position position
Definition: ComputeQM.C:66
void NAMD_bug(const char *err_msg)
Definition: common.C:129
void recvFullQM(QMCoordMsg *)
Definition: ComputeQM.C:1281
#define QMCHRGMULLIKEN
Definition: Molecule.h:134
int get_qmPCFreq()
Definition: Molecule.h:862
#define QMPCTYPE_IGNORE
Definition: ComputeQM.C:197
gridSize z
int Bool
Definition: common.h:133
int numRealPntChrgs
Definition: ComputeQM.C:325
Force * f[maxNumForces]
Definition: PatchTypes.h:67
int write_dcdstep(int fd, int N, float *X, float *Y, float *Z, double *cell)
Definition: dcdlib.C:736
Real * get_qmAtmChrg()
Definition: Molecule.h:823
void clear()
Definition: ResizeArray.h:87
const int *const * get_qmMMBond()
Definition: Molecule.h:837
bool operator==(const ComputeQMPntChrg &ref)
Definition: ComputeQM.C:162
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
int get_qmLSSFreq()
Definition: Molecule.h:848
int sourceNode
Definition: ComputeQM.C:105
float charge
Definition: ComputeQM.h:105
QMForce * force
Definition: ComputeQM.C:191
char prepProc[256]
Definition: ComputeQM.C:338
BigReal x
Definition: Vector.h:66
idIndxStr(int newID, int newIndx)
Definition: ComputeQM.C:212
std::vector< ComputeQMPntChrg > QMPCVec
Definition: ComputeQM.C:398
bool custom_ComputeQMAtom_Less(const ComputeQMAtom a, const ComputeQMAtom b)
Definition: ComputeQM.C:94
Real qmGrpID
Definition: ComputeQM.C:69
int index(const Elem &elem)
Definition: SortedArray.h:75
BigReal PMEEwaldCoefficient
bool operator<(const ComputeQMPntChrg &ref)
Definition: ComputeQM.C:159
const Real * get_cSMDcoffs()
Definition: Molecule.h:875
int numAtoms
Definition: Molecule.h:557
void NAMD_die(const char *err_msg)
Definition: common.C:85
char element[3]
Definition: ComputeQM.C:301
Force force
Definition: NamdTypes.h:202
static char * FORMAT(BigReal X)
Definition: ComputeQM.C:367
ComputeQMPntChrg * coord
Definition: ComputeQM.C:171
#define QMLENTYPE
Definition: ComputeQM.C:55
ConfigList * configList
Definition: Node.h:179
void calcORCA(QMGrpCalcMsg *)
Definition: ComputeQM.C:3649
LSSSubsDat * lssDat
Definition: ComputeQM.C:192
Force force
Definition: ComputeQM.h:103
const int *const * get_qmMMChargeTarget()
Definition: Molecule.h:842
int * get_qmCustPCSizes()
Definition: Molecule.h:865
int write_dcdheader(int fd, const char *filename, int N, int NFILE, int NPRIV, int NSAVC, int NSTEP, double DELTA, int with_unitcell)
Definition: dcdlib.C:915
void processFullQM(QMCoordMsg *)
Definition: ComputeQM.C:1292
ComputeQM(ComputeID c)
Definition: ComputeQM.C:599
const char *const * get_qmElements()
Definition: Molecule.h:827
#define QMATOMTYPE_QM
Definition: ComputeQM.C:196
int add(const Elem &elem)
Definition: ResizeArray.h:97
Bool qmMOPACAddConfigChrg
int pcSelMode
Definition: ComputeQM.C:109
void calcUSR(QMGrpCalcMsg *)
Definition: ComputeQM.C:4511
void saveResults(QMForceMsg *)
Definition: ComputeQM.C:2673
BlockRadixSort::TempStorage sort
int resIndx
Definition: ComputeQM.C:359
BigReal length2(void) const
Definition: Vector.h:173
void calcMOPAC(QMGrpCalcMsg *)
Definition: ComputeQM.C:2801
const Real * get_cSMDKs()
Definition: Molecule.h:873
void recvForce(QMForceMsg *)
Definition: ComputeQM.C:2662
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ cudaTextureObject_t const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cutoff2
#define simParams
Definition: Output.C:127
short vdwType
Definition: NamdTypes.h:55
StringList * next
Definition: ConfigList.h:49
void resize(int i)
Definition: ResizeArray.h:84
virtual void initialize()
char secProc[256]
Definition: ComputeQM.C:338
char * data
Definition: ConfigList.h:48
BigReal y
Definition: Vector.h:66
#define QMPCSCALESWITCH
Definition: ComputeQM.C:63
bool operator==(const lssDistSort &ref)
Definition: ComputeQM.C:274
BigReal virial[3][3]
Definition: ComputeQM.C:188
void initialize()
Definition: ComputeQM.C:615
Mass mass
Definition: NamdTypes.h:108
StringList * find(const char *name) const
Definition: ConfigList.C:341
std::map< Real, std::pair< Position, BigReal > > GrpDistMap
Definition: ComputeQM.C:1289
std::pair< int, LSSDataStr > atmLSSData
Definition: ComputeQM.C:392
Elem * find(const Elem &elem)
#define TIMEFACTOR
Definition: common.h:48
int numAtoms
Definition: ComputeQM.C:106
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
int * get_qmCustomPCIdxs()
Definition: Molecule.h:866
gridSize y
int homeIndx
Definition: ComputeQM.h:104
const Real * get_cSMDVels()
Definition: Molecule.h:874
int * get_qmLSSSize()
Definition: Molecule.h:845
ComputeQMAtom(const ComputeQMAtom &ref)
Definition: ComputeQM.C:83
BigReal dielectric
void submit(void)
Definition: ReductionMgr.h:323
int size(void) const
Definition: ResizeArray.h:127
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:83
int numAllPntChrgs
Definition: ComputeQM.C:326
int bondIndx
Definition: ComputeQM.C:350
#define QMPCSCHEMEZERO
Definition: ComputeQM.C:60
#define QMPCTYPE_EXTRA
Definition: ComputeQM.C:199
#define PCMODECUSTOMSEL
Definition: ComputeQM.C:101
#define QMLSSMODEDIST
Definition: Molecule.h:126
gridSize x
BigReal virial[3][3]
Definition: ComputeQM.C:179
void close_dcd_write(int fd)
Definition: dcdlib.C:1063
LSSDataStr(int newResIndx, Mass newMass)
Definition: ComputeQM.C:361
int const *const * get_cSMDindex()
Definition: Molecule.h:871
Mass * get_qmLSSRefMass()
Definition: Molecule.h:852
char * configLines
Definition: ComputeQM.C:340
Position pos
Definition: ComputeQM.C:345
char qmBaseDir[NAMD_FILENAME_BUFFER_SIZE]
int load(const Elem &elem)
Definition: SortedArray.h:50
ComputeQMPntChrg(const ComputeQMPntChrg &ref)
Definition: ComputeQM.C:148
Molecule * molecule
Definition: Node.h:176
const BigReal * get_qmDummyBondVal()
Definition: Molecule.h:835
bool secProcOn
Definition: ComputeQM.C:329
#define QMLSSQMRES
Definition: ComputeQM.C:201
Bool get_qmReplaceAll()
Definition: Molecule.h:855
int numPCIndxs
Definition: ComputeQM.C:108
ResizeArrayIter< T > begin(void) const
Vector unit(void) const
Definition: Vector.h:182
QMAtomData(Position posInit, float chrgInit, int idInit, int bountToIndxInit, int newType, char *elementInit, Real newDist)
Definition: ComputeQM.C:305
#define QMLSSCLASSICALRES
Definition: ComputeQM.C:202
std::map< int, LSSDataStr > LSSDataMap
Definition: ComputeQM.C:393
double BigReal
Definition: common.h:114
int get_qmLSSResSize()
Definition: Molecule.h:849
const int * get_qmGrpSizes()
Definition: Molecule.h:829
std::vector< ComputeQMAtom > QMAtmVec
Definition: ComputeQM.C:399
#define QMCHRGNONE
Definition: Molecule.h:133
BigReal energyCorr
Definition: ComputeQM.C:178
for(int i=0;i< n1;++i)
bool PMEOn
Definition: ComputeQM.C:186