ComputeQM.C

Go to the documentation of this file.
00001 
00007 #include "dcdlib.h"
00008 
00009 // #define DEBUGM
00010 // #define DEBUG_QM
00011 
00012 #ifdef DEBUG_QM
00013   #define DEBUGM
00014 #endif
00015 
00016 #define MIN_DEBUG_LEVEL 1
00017 #include "Debug.h"
00018 
00019 #include "InfoStream.h"
00020 #include "Node.h"
00021 #include "PatchMap.h"
00022 #include "PatchMap.inl"
00023 #include "AtomMap.h"
00024 #include "ComputeQM.h"
00025 #include "ComputeQMMgr.decl.h"
00026 #include "PatchMgr.h"
00027 #include "Molecule.h"
00028 #include "ReductionMgr.h"
00029 #include "ComputeMgr.h"
00030 #include "ComputeMgr.decl.h"
00031 #include "SimParameters.h"
00032 #include "WorkDistrib.h"
00033 #include "varsizemsg.h"
00034 #include <stdlib.h>
00035 #include <stdio.h>
00036 #include <errno.h>
00037 #include <algorithm>
00038 
00039 #include "ComputePme.h"
00040 #include "ComputePmeMgr.decl.h"
00041 
00042 #include <fstream>
00043 #include <iomanip>
00044 
00045 #if defined(WIN32) && !defined(__CYGWIN__)
00046 #include <direct.h>
00047 #define mkdir(X,Y) _mkdir(X)
00048 #define S_ISDIR(X) ((X) & S_IFDIR)
00049 #endif
00050 
00051 #ifndef SQRT_PI
00052 #define SQRT_PI 1.7724538509055160273 /* mathematica 15 digits*/
00053 #endif
00054 
00055 #define QMLENTYPE 1
00056 #define QMRATIOTYPE 2
00057 
00058 #define QMPCSCHEMENONE 1
00059 #define QMPCSCHEMEROUND 2
00060 #define QMPCSCHEMEZERO 3
00061 
00062 #define QMPCSCALESHIFT 1
00063 #define QMPCSCALESWITCH 2
00064 
00065 struct ComputeQMAtom {
00066     Position position;
00067     float charge;
00068     int id;
00069     Real qmGrpID;
00070     int homeIndx;
00071     int vdwType;
00072     ComputeQMAtom() : position(0), charge(0), id(-1), qmGrpID(-1), 
00073                         homeIndx(-1), vdwType(0) {}
00074     ComputeQMAtom(Position posInit, float chrgInit, int idInit, 
00075                   Real qmInit, int hiInit, int newvdwType) {
00076         position = posInit;
00077         charge = chrgInit;
00078         id = idInit;
00079         qmGrpID = qmInit;
00080         homeIndx = hiInit;
00081         vdwType = newvdwType;
00082     }
00083     ComputeQMAtom(const ComputeQMAtom &ref) {
00084         position = ref.position;
00085         charge = ref.charge;
00086         id = ref.id;
00087         qmGrpID = ref.qmGrpID;
00088         homeIndx = ref.homeIndx;
00089         vdwType = ref.vdwType;
00090     }
00091 };
00092 
00093 // sort using a custom function
00094 bool custom_ComputeQMAtom_Less(const ComputeQMAtom a,const ComputeQMAtom b)
00095 {   
00096     return a.id < b.id;
00097 }
00098 
00099 #define PCMODEUPDATESEL 1
00100 #define PCMODEUPDATEPOS 2
00101 #define PCMODECUSTOMSEL 3
00102 
00103 class QMCoordMsg : public CMessage_QMCoordMsg {
00104 public:
00105   int sourceNode;
00106   int numAtoms;
00107   int timestep;
00108   int numPCIndxs;
00109   int pcSelMode;
00110   ComputeQMAtom *coord;
00111   int *pcIndxs;
00112 };
00113 
00114 struct pntChrgDist {
00115     int index;
00116     Real dist;
00117     pntChrgDist() : index(-1), dist(0) {};
00118     pntChrgDist(int newIndx, Real newDist) {
00119         index = newIndx;
00120         dist = newDist;
00121     }
00122     int operator <(const pntChrgDist& v) {return dist < v.dist;}
00123 };
00124 
00125 struct ComputeQMPntChrg {
00126     Position position;
00127     float charge;
00128     int id;
00129     Real qmGrpID;
00130     int homeIndx;
00131     Real dist;
00132     Mass mass;
00133     int vdwType;
00134     ComputeQMPntChrg() : position(0), charge(0), id(-1), qmGrpID(-1), 
00135         homeIndx(-1), dist(0), mass(0), vdwType(0) {}
00136     ComputeQMPntChrg(Position posInit, float chrgInit, int idInit, 
00137                      Real qmInit, int hiInit, Real newDist, 
00138                      Mass newM, int newvdwType) {
00139         position = posInit;
00140         charge = chrgInit;
00141         id = idInit;
00142         qmGrpID = qmInit;
00143         homeIndx = hiInit;
00144         dist = newDist;
00145         mass = newM;
00146         vdwType = newvdwType;
00147     }
00148     ComputeQMPntChrg(const ComputeQMPntChrg &ref) {
00149         position = ref.position;
00150         charge = ref.charge;
00151         id = ref.id;
00152         qmGrpID = ref.qmGrpID;
00153         homeIndx = ref.homeIndx;
00154         dist = ref.dist;
00155         mass = ref.mass;
00156         vdwType = ref.vdwType;
00157     }
00158     
00159     bool operator <(const ComputeQMPntChrg& ref) {
00160         return (id < ref.id);
00161     }
00162     bool operator ==(const ComputeQMPntChrg& ref) {
00163         return (id == ref.id) ;
00164     }
00165 };
00166 
00167 class QMPntChrgMsg : public CMessage_QMPntChrgMsg {
00168 public:
00169   int sourceNode;
00170   int numAtoms;
00171   ComputeQMPntChrg *coord;
00172 };
00173 
00174 class QMGrpResMsg : public CMessage_QMGrpResMsg {
00175 public:
00176   int grpIndx;
00177   BigReal energyOrig; // Original QM Energy
00178   BigReal energyCorr; // Corrected energy due to PME
00179   BigReal virial[3][3];
00180   int numForces;
00181   QMForce *force;
00182 };
00183 
00184 class QMForceMsg : public CMessage_QMForceMsg {
00185 public:
00186   bool PMEOn;
00187   BigReal energy;
00188   BigReal virial[3][3];
00189   int numForces;
00190   int numLssDat;
00191   QMForce *force;
00192   LSSSubsDat *lssDat;
00193 };
00194 
00195 # define QMATOMTYPE_DUMMY 0
00196 # define QMATOMTYPE_QM 1
00197 # define QMPCTYPE_IGNORE 0
00198 # define QMPCTYPE_CLASSICAL 1
00199 # define QMPCTYPE_EXTRA 2
00200 
00201 #define QMLSSQMRES 1
00202 #define QMLSSCLASSICALRES 2
00203 
00204 struct idIndxStr {
00205     int ID; // Global atom ID
00206     int indx; // Atom index in the vector received from other ranks.
00207     
00208     idIndxStr() {
00209         ID = -1;
00210         indx = -1;
00211     };
00212     idIndxStr(int newID, int newIndx) {
00213         ID = newID;
00214         indx = newIndx;
00215     }
00216     idIndxStr(const idIndxStr& ref) {
00217         ID = ref.ID ;
00218         indx = ref.indx;
00219     }
00220     
00221     idIndxStr &operator=(const idIndxStr& ref) {
00222         ID = ref.ID ;
00223         indx = ref.indx;
00224         return *this;
00225     }
00226     
00227     bool operator==(const idIndxStr& ref) const {
00228         return (ID == ref.ID && indx == ref.indx);
00229     }
00230     bool operator!=(const idIndxStr& ref) const {
00231         return !(*this == ref);
00232     }
00233     
00234     // We sort the atoms based on their global ID so that the same
00235     // order is kept when replacing classical and QM residues.
00236     // We are assuming here that every time a residue is added to a PDB
00237     // file, its atoms are placed in the same order.
00238     bool operator<(const idIndxStr& ref) {return ID < ref.ID;}
00239 };
00240 
00241 struct lssDistSort {
00242     int type;
00243     Real dist;
00244 //     std::vector<int> idVec;
00245 //     std::vector<int> indxVec;
00246     SortedArray<idIndxStr> idIndx;
00247     
00248     lssDistSort() {
00249         idIndx.resize(0);
00250     };
00251     lssDistSort(int newType, Real newDist) {
00252         type = newType;
00253         dist = newDist;
00254         idIndx.resize(0);
00255     }
00256     lssDistSort(const lssDistSort& ref) {
00257         type = ref.type;
00258         dist = ref.dist;
00259         idIndx.resize(idIndx.size());
00260         for (int i=0; i<ref.idIndx.size(); i++) idIndx.insert(ref.idIndx[i]) ;
00261         idIndx.sort();
00262     }
00263     
00264     lssDistSort& operator=(const lssDistSort& ref) {
00265         type = ref.type;
00266         dist = ref.dist;
00267         idIndx.resize(idIndx.size());
00268         for (int i=0; i<ref.idIndx.size(); i++) idIndx.insert(ref.idIndx[i]) ;
00269         idIndx.sort();
00270         
00271         return *this;
00272     }
00273     
00274     bool operator==(const lssDistSort& ref) {
00275         bool returnVal = true;
00276         
00277         if (! (type == ref.type && dist == ref.dist))
00278             return false;
00279         
00280         if (idIndx.size() != ref.idIndx.size())
00281             return false;
00282         
00283         for (int i=0; i<ref.idIndx.size(); i++) {
00284             if (idIndx[i] != ref.idIndx[i])
00285                 return false;
00286         }
00287         
00288         return returnVal;
00289     }
00290     
00291     bool operator<(const lssDistSort& ref) {return dist < ref.dist;}
00292     
00293 } ;
00294 
00295 struct QMAtomData {
00296     Position position;
00297     float charge;
00298     int id;
00299     int bountToIndx; // So as to implement an "exclude 1-2" functionality.
00300     int type;
00301     char element[3];
00302     Real dist;
00303     QMAtomData() : position(0), charge(0), id(-1), bountToIndx(-1), 
00304                    type(-1), dist(0) {}
00305     QMAtomData(Position posInit, float chrgInit, int idInit, 
00306                int bountToIndxInit, int newType, 
00307                char *elementInit, Real newDist) {
00308         position = posInit;
00309         charge = chrgInit;
00310         id = idInit;
00311         bountToIndx = bountToIndxInit;
00312         type = newType;
00313         strncpy(element,elementInit,3);
00314         dist = newDist;
00315     }
00316 };
00317 
00318 class QMGrpCalcMsg: public CMessage_QMGrpCalcMsg {
00319 public:
00320     int timestep;
00321     int grpIndx;
00322     int peIter;
00323     int numQMAtoms ;
00324     int numAllAtoms ; // Including dummy atoms.
00325     int numRealPntChrgs;
00326     int numAllPntChrgs; // Inlcuding point charges created to handle QM-MM bonds.
00327     Real charge, multiplicity;
00328     BigReal constants;
00329     bool secProcOn ;
00330     bool prepProcOn ;
00331     bool PMEOn;
00332     bool switching ;
00333     int switchType;
00334     Real cutoff, swdist;
00335     int pcScheme ;
00336     BigReal PMEEwaldCoefficient;
00337     int qmAtmChrgMode;
00338     char baseDir[256], execPath[256], secProc[256], prepProc[256];
00339     QMAtomData *data;
00340     char *configLines;
00341 };
00342 
00343 struct dummyData {
00344     // Position of dummy atom.
00345     Position pos ;
00346     // Indicates to which QM atom is this dummy atom bound.
00347     // This is indexed in the context of the number of QM atoms in a QM group.
00348     int qmInd ;
00349     // Indicates the index of this bond, used to get the element of the dummy atom.
00350     int bondIndx;
00351     dummyData(Position newPos, int newQmInd, int newIndx) { 
00352         pos = newPos;
00353         qmInd = newQmInd;
00354         bondIndx = newIndx;
00355     }
00356 } ;
00357 
00358 struct LSSDataStr {
00359     int resIndx ;
00360     Mass mass ;
00361     LSSDataStr(int newResIndx, Mass newMass) {
00362         resIndx = newResIndx;
00363         mass = newMass;
00364     }
00365 } ;
00366 
00367 static char *FORMAT(BigReal X)
00368 {
00369   static char tmp_string[25];
00370   const double maxnum = 9999999999.9999;
00371   if ( X > maxnum ) X = maxnum;
00372   if ( X < -maxnum ) X = -maxnum;
00373   sprintf(tmp_string," %14.4f",X); 
00374   return tmp_string;
00375 }
00376 
00377 static char *FORMAT(const char *X)
00378 {
00379   static char tmp_string[25];
00380   sprintf(tmp_string," %14s",X); 
00381   return tmp_string;
00382 }
00383 
00384 static char *QMETITLE(int X)
00385 {
00386   static char tmp_string[21];
00387   sprintf(tmp_string,"QMENERGY: %7d",X); 
00388   return  tmp_string;
00389 }
00390 
00391 
00392 typedef std::pair<int, LSSDataStr> atmLSSData;
00393 typedef std::map<int, LSSDataStr> LSSDataMap;
00394 
00395 typedef std::pair<int, Mass> refLSSData;
00396 typedef std::map<int, Mass> LSSRefMap;
00397 
00398 typedef std::vector<ComputeQMPntChrg> QMPCVec ;
00399 typedef std::vector<ComputeQMAtom> QMAtmVec ;
00400 
00401 class ComputeQMMgr : public CBase_ComputeQMMgr {
00402 public:
00403     ComputeQMMgr();
00404     ~ComputeQMMgr();
00405 
00406     void setCompute(ComputeQM *c) { QMCompute = c; }
00407     
00408     // Function called on pe ZERO by all pe's.
00409     // Receives partial QM coordinates form the QM atoms in each pe and aggregates
00410     // them all in a single structure.
00411     // This function initializes some variables and allcoates memory.
00412     void recvPartQM(QMCoordMsg*) ;
00413     
00414     void recvFullQM(QMCoordMsg*) ;
00415     
00416     void recvPntChrg(QMPntChrgMsg*) ;
00417     
00418     void calcMOPAC(QMGrpCalcMsg *) ;
00419     void calcORCA(QMGrpCalcMsg *) ;
00420     void calcUSR(QMGrpCalcMsg *) ;
00421     
00422     void storeQMRes(QMGrpResMsg *) ;
00423     
00424     void procQMRes();
00425     
00426     void recvForce(QMForceMsg*) ;
00427     
00428     SortedArray<LSSSubsDat> &get_subsArray() { return subsArray; } ;
00429 private:
00430     ComputeQMMgr_SDAG_CODE
00431     
00432     CProxy_ComputeQMMgr QMProxy;
00433     ComputeQM *QMCompute;
00434 
00435     int numSources;
00436     int numArrivedQMMsg, numArrivedPntChrgMsg, numRecQMRes;
00437     QMCoordMsg **qmCoordMsgs;
00438     QMPntChrgMsg **pntChrgCoordMsgs;
00439     
00440     std::vector<int> qmPEs ;
00441     
00442     ComputeQMAtom *qmCoord;
00443     QMPCVec pntChrgMsgVec;
00444     QMForce *force;
00445 
00446     int numQMGrps;
00447 
00448     int numAtoms;
00449 
00450     int qmAtmIter;
00451 
00452     int numQMAtms;
00453     
00454     Bool noPC ;
00455     int meNumMMIndx;
00456     
00457     int qmPCFreq;
00458     // In case we are skiping steps between point charge re-assignment (i.e. qmPCFreq > 0),
00459     // we keep a list in rank zero whihch is replaced every X steps, and 
00460     // then re-sent to all PEs every "X+1" steps. This is done because atoms
00461     // may move between pathes and PEs between PC reassignments, if they are 
00462     // sufficiently long or unlucky.
00463     int *pcIDList ;
00464     int pcIDListSize;
00465     Bool resendPCList;
00466     SortedArray<ComputeQMPntChrg> pcDataSort;
00467     
00468     int replaceForces;
00469     int bondValType;
00470     SimParameters * simParams;
00471     Molecule *molPtr;
00472     // ID for each QM group.
00473     Real *grpID;
00474     Real *grpChrg, *grpMult;
00475     
00476     Bool qmLSSOn ;
00477     int qmLSSFreq;
00478     int qmLSSResSize;
00479     int *qmLSSSize;
00480     int *qmLSSIdxs;
00481     Mass *qmLSSMass;
00482     int lssTotRes;
00483     // This will hold one map per QM group. Each will map an atom ID with
00484     // the local residue index for the solvent molecule it belongs to, 
00485     // and its atomic mass (in case COM calculation is being performed).
00486     ResizeArray<LSSDataMap> grpIDResNum ;
00487     int *qmLSSRefIDs;
00488     int *qmLSSRefSize;
00489     // This will hold one map per QM group. Each will map an atom ID with
00490     // its atomic mass. This is used to calculate the reference COM to which
00491     // solvent molecules will be compared for distance calculation.
00492     ResizeArray<LSSRefMap> lssGrpRefMass ;
00493     // Current positions of each LSS RESIDUE
00494     Position *lssPos;
00495     Mass lssResMass;
00496     SortedArray<LSSSubsDat> subsArray;
00497     
00498     BigReal ewaldcof, pi_ewaldcof;
00499     
00500     // Accumulates the enery of different QM groups. This is the energy passed to NAMD.
00501     BigReal totalEnergy;
00502     BigReal totVirial[3][3];
00503 
00504     int dcdOutFile, dcdPosOutFile;
00505     Real *outputData ;
00506     int timeStep ;
00507     
00508     void procBonds(int numBonds,
00509                      const int *const qmGrpBondsGrp,
00510                      const int *const qmMMBondedIndxGrp,
00511                      const int *const *const chargeTarget,
00512                      const int *const numTargs,
00513                      const QMPCVec grpPntChrgVec,
00514                      QMPCVec &grpAppldChrgVec) ;
00515     
00516     void pntChrgSwitching(QMGrpCalcMsg* msg, QMAtomData *pcPpme) ;
00517     
00518     void lssPrepare() ;
00519     void lssUpdate(int grpIter, QMAtmVec &grpQMAtmVec, QMPCVec &grpPntChrgVec);
00520         
00521         void calcCSMD(int grpIndx,int numQMAtoms, const QMAtomData *atmP, Force *resForce) ;
00522         Bool cSMDon;
00523 
00524     int cSMDnumInstances, cSMDInitGuides;
00525     // Index of Conditional SMD guides per QM group
00526     int const * const * cSMDindex;
00527     // Num instances per QM group
00528     int const * cSMDindxLen;
00529     // Positions of Conditional SMD guides
00530     Position* cSMDguides;
00531     // Atom indices for Origin and Target of cSMD
00532     int const * const * cSMDpairs;
00533     // Spring constants for cSMD
00534     Real const * cSMDKs;
00535     // Speed of movement of guide particles for cSMD.
00536     Real const * cSMDVels;
00537     // Distance cutoff for guide particles for cSMD.
00538     Real const * cSMDcoffs;
00539     // Vector where we store all calculated cSMD forces
00540     Force * cSMDForces ;
00541     
00542     #ifdef DEBUG_QM
00543     void Write_PDB(std::string Filename, const QMGrpCalcMsg *dataMsg);
00544     void Write_PDB(std::string Filename, const QMCoordMsg *dataMsg);
00545     #endif
00546 };
00547 
00548 ComputeQMMgr::ComputeQMMgr() :
00549   QMProxy(thisgroup), QMCompute(0), numSources(0), numArrivedQMMsg(0), 
00550   numArrivedPntChrgMsg(0), numRecQMRes(0), qmCoordMsgs(0), pntChrgCoordMsgs(0),
00551   qmCoord(0), force(0), numAtoms(0), dcdOutFile(0), outputData(0), pcIDListSize(0) {
00552       
00553   CkpvAccess(BOCclass_group).computeQMMgr = thisgroup;
00554   
00555 }
00556 
00557 ComputeQMMgr::~ComputeQMMgr() {
00558     
00559     if (qmCoordMsgs != 0) {
00560         for ( int i=0; i<numSources; ++i ) {
00561             if (qmCoordMsgs[i] != 0)
00562                 delete qmCoordMsgs[i];
00563         }
00564         delete [] qmCoordMsgs;
00565     }
00566     
00567     if (pntChrgCoordMsgs != 0) {
00568         for ( int i=0; i<numSources; ++i ) {
00569             if (pntChrgCoordMsgs[i] != 0)
00570                 delete pntChrgCoordMsgs[i];
00571         }
00572         delete [] pntChrgCoordMsgs;
00573     }
00574     
00575     if (qmCoord != 0)
00576         delete [] qmCoord;
00577     
00578     if (force != 0)
00579         delete [] force;
00580     
00581     
00582     if (dcdOutFile != 0)
00583         close_dcd_write(dcdOutFile);
00584     
00585     if (dcdPosOutFile != 0)
00586         close_dcd_write(dcdPosOutFile);
00587     
00588     if (outputData != 0)
00589         delete [] outputData;
00590     
00591     if (lssPos != NULL)
00592         delete [] lssPos;
00593 }
00594 
00595 SortedArray<LSSSubsDat> &lssSubs(ComputeQMMgr *mgr) { 
00596     return mgr->get_subsArray(); 
00597 } ;
00598 
00599 ComputeQM::ComputeQM(ComputeID c) :
00600   ComputeHomePatches(c), oldForces(0)
00601 {
00602   CProxy_ComputeQMMgr::ckLocalBranch(
00603         CkpvAccess(BOCclass_group).computeQMMgr)->setCompute(this);
00604 
00605   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00606   
00607 }
00608 
00609 ComputeQM::~ComputeQM()
00610 {
00611     if (oldForces != 0)
00612         delete [] oldForces;
00613 }
00614 
00615 void ComputeQM::initialize()
00616 {
00617     ComputeHomePatches::initialize();
00618     
00619     // Get some pointers that will be used in the future,
00620     simParams = Node::Object()->simParameters ;
00621     molPtr = Node::Object()->molecule;
00622     // Future comes fast...
00623     qmAtomGroup = molPtr->get_qmAtomGroup() ;
00624     
00625     noPC = molPtr->get_noPC();
00626     if (noPC) {
00627         meNumMMIndx = molPtr->get_qmMeNumBonds();
00628         meMMindx = molPtr->get_qmMeMMindx();
00629         meQMGrp =molPtr->get_qmMeQMGrp();
00630         
00631         for (int i=0; i<meNumMMIndx; i++)
00632             meQMBonds.add(meMMQMGrp(meMMindx[i], meQMGrp[i]));
00633     }
00634     
00635     numQMAtms = molPtr->get_numQMAtoms();
00636     qmAtmChrg = molPtr->get_qmAtmChrg() ;
00637     qmAtmIndx = molPtr->get_qmAtmIndx() ;
00638     
00639     numQMGrps = molPtr->get_qmNumGrps();
00640     qmGrpIDArray = molPtr->get_qmGrpID() ;
00641     
00642     cutoff = simParams->cutoff;
00643     
00644     customPC = simParams->qmCustomPCSel;
00645     if (customPC) {
00646         
00647         customPCLists.resize(numQMGrps);
00648         
00649         int totCustPCs = molPtr->get_qmTotCustPCs();
00650         const int *custPCSizes = molPtr->get_qmCustPCSizes();
00651         const int *customPCIdxs = molPtr->get_qmCustomPCIdxs();
00652         
00653         int minI = 0, maxI = 0, grpIter = 0;
00654         
00655         while (grpIter < numQMGrps) {
00656             
00657             maxI += custPCSizes[grpIter];
00658             
00659             for (size_t i=minI; i < maxI; i++) {
00660                 
00661                 customPCLists[grpIter].add( customPCIdxs[i] ) ;
00662             }
00663             
00664             // Minimum index gets the size of the current group
00665             minI += custPCSizes[grpIter];
00666             // Group iterator gets incremented
00667             grpIter++;
00668             
00669         }
00670     }
00671     
00672 }
00673 
00674 
00675 void ComputeQM::doWork()
00676 {
00677     CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
00678     
00679     ResizeArrayIter<PatchElem> ap(patchList);
00680     
00681         int timeStep;
00682         
00683     #ifdef DEBUG_QM
00684     DebugM(4,"----> Initiating QM work on rank " << CkMyPe() <<
00685     " with " << patchList.size() << " patches." << std::endl );
00686     #endif
00687     
00688     patchData.clear();
00689     
00690     // Determines hou many QM atoms are in this node.
00691     int numLocalQMAtoms = 0, localMM1atoms = 0;
00692     for (ap = ap.begin(); ap != ap.end(); ap++) {
00693         CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
00694         int localNumAtoms = (*ap).p->getNumAtoms();
00695         
00696         patchData.push_back(patchDataStrc((*ap).positionBox,
00697                     (*ap).positionBox->open(),
00698                     (*ap).p) );
00699         
00700         for(int i=0; i<localNumAtoms; ++i) {
00701             if ( qmAtomGroup[xExt[i].id] > 0 ) {
00702                 numLocalQMAtoms++;
00703             }
00704             else if (meNumMMIndx > 0) {
00705                 // If we have a mechanical embedding simulation with no 
00706                 // point charges, but with QM-MM bonds, we look for the MM1 atoms 
00707                 // here and pass them directly. No need for an extra round of 
00708                 // comms just to get atoms we already know we need.
00709                 
00710                 auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
00711                 
00712                 if (retIt != NULL) {
00713                     localMM1atoms++;
00714                 }
00715             }
00716         }
00717         
00718         timeStep = (*ap).p->flags.step ;
00719     }
00720     
00721     DebugM(4, "Node " << CkMyPe() << " has " << numLocalQMAtoms
00722     << " QM atoms." << std::endl) ;
00723     #ifdef DEBUG_QM
00724     if (localMM1atoms > 0)
00725         DebugM(4, "Node " << CkMyPe() << " has " << localMM1atoms
00726             << " MM1 atoms." << std::endl) ;
00727     #endif
00728     // Send QM atoms
00729     
00730     // This pointer will be freed in rank zero.
00731     QMCoordMsg *msg = new (numLocalQMAtoms + localMM1atoms,0, 0) QMCoordMsg;
00732     msg->sourceNode = CkMyPe();
00733     msg->numAtoms = numLocalQMAtoms + localMM1atoms;
00734     msg->timestep = timeStep;
00735     ComputeQMAtom *data_ptr = msg->coord;
00736     
00737     // Build a message with coordinates, charge, atom ID and QM group
00738     // for all QM atoms in the node, then send it to node zero.
00739     int homeCounter = 0;
00740     for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
00741         
00742         CompAtom *x = (*pdIt).compAtomP;
00743         const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
00744         int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
00745         const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
00746         
00747         for(int i=0; i<localNumAtoms; ++i) {
00748             
00749             if ( qmAtomGroup[xExt[i].id] > 0 ) {
00750                 
00751                 Real charge = 0;
00752                 
00753                 for (int qi=0; qi<numQMAtms; qi++) {
00754                     if (qmAtmIndx[qi] == xExt[i].id) {
00755                         charge = qmAtmChrg[qi];
00756                         break;
00757                     }
00758                 }
00759                 
00760                 data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position, 
00761                                                           fullAtms[i].transform) ;
00762 //                 data_ptr->charge = x[i].charge;
00763                 data_ptr->charge = charge;
00764                 data_ptr->id = xExt[i].id;
00765                 data_ptr->qmGrpID = qmAtomGroup[xExt[i].id] ;
00766                 data_ptr->homeIndx = homeCounter;
00767                 data_ptr->vdwType = fullAtms[i].vdwType;
00768                 
00769                 ++data_ptr;
00770             }
00771             else if (meNumMMIndx > 0) {
00772                 // If we have a mechanical embedding simulation with no 
00773                 // point charges, but with QM-MM bonds, we look for the MM1 atoms 
00774                 // connected here and pass them directly.
00775                 
00776                 auto retIt = meQMBonds.find(meMMQMGrp(xExt[i].id));
00777                 
00778                 if (retIt != NULL) {
00779                     
00780                     DebugM(3,"Found atom " << retIt->mmIndx << "," << retIt->qmGrp << "\n" );
00781                     
00782                     data_ptr->position = patchList[0].p->lattice.reverse_transform(x[i].position, 
00783                                                           fullAtms[i].transform) ;
00784                     // We force the charge to be zero since we would only get
00785                     // here if mechanical embeding was selected.
00786                     data_ptr->charge = 0;
00787                     data_ptr->id = xExt[i].id;
00788                     data_ptr->qmGrpID = retIt->qmGrp ;
00789                     data_ptr->homeIndx = homeCounter;
00790                     
00791                     // We re-use this varaible to indicate this is an MM1 atom.
00792                     data_ptr->vdwType = -1;
00793                     
00794                     ++data_ptr;
00795                     
00796                 }
00797                 
00798             }
00799                 
00800             homeCounter++;
00801             
00802         }
00803         
00804     }
00805     
00806     if (noPC) {
00807         for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
00808             CompAtom *x = (*pdIt).compAtomP;
00809             (*pdIt).posBoxP->close(&x);
00810         }
00811     }
00812     
00813     QMProxy[0].recvPartQM(msg);
00814     
00815 }
00816 
00817 
00818 void ComputeQMMgr::recvPartQM(QMCoordMsg*msg)
00819 {
00820     // In the first (ever) step of the simulation, we allocate arrays that
00821     // will be used throughout the simulation, and get some important numbers.
00822     if ( ! numSources ) {
00823         DebugM(4,"Initializing ComputeQMMgr variables." << std::endl);
00824         numSources = (PatchMap::Object())->numNodesWithPatches();
00825         
00826         DebugM(4,"There are " << numSources << " nodes with patches." << std::endl);
00827         qmCoordMsgs = new QMCoordMsg*[numSources];
00828         for ( int i=0; i<numSources; ++i ) { 
00829             qmCoordMsgs[i] = 0;
00830         }
00831         
00832         // Prepares the allocation for the recvPntChrg function.
00833         
00834         DebugM(4,"Getting data from molecule and simParameters." << std::endl);
00835         
00836         molPtr = Node::Object()->molecule;
00837         simParams = Node::Object()->simParameters;
00838         
00839         numAtoms = molPtr->numAtoms;
00840         force = new QMForce[numAtoms];
00841         
00842         numQMAtms = molPtr->get_numQMAtoms();
00843         qmAtmIter = 0;
00844         
00845         noPC = simParams->qmNoPC ;
00846         meNumMMIndx = molPtr->get_qmMeNumBonds();
00847         if (noPC && meNumMMIndx == 0) {
00848             pntChrgCoordMsgs = NULL;
00849         }
00850         else {
00851             pntChrgCoordMsgs = new QMPntChrgMsg*[numSources];
00852             for ( int i=0; i<numSources; ++i ) { 
00853                 pntChrgCoordMsgs[i] = 0;
00854             }
00855         }
00856         
00857         qmPCFreq = molPtr->get_qmPCFreq();
00858         resendPCList = false;
00859         
00860         grpID = molPtr->get_qmGrpID() ;
00861         bondValType = simParams->qmBondValType;
00862         
00863         numQMGrps = molPtr->get_qmNumGrps();
00864         
00865         grpChrg = molPtr->get_qmGrpChrg() ;
00866         
00867         grpMult = molPtr->get_qmGrpMult() ;
00868         
00869         qmLSSOn = simParams->qmLSSOn ;
00870         if (qmLSSOn) {
00871             qmLSSFreq = molPtr->get_qmLSSFreq() ;
00872             qmLSSSize = molPtr->get_qmLSSSize() ;
00873             qmLSSIdxs = molPtr->get_qmLSSIdxs() ;
00874             qmLSSMass = molPtr->get_qmLSSMass() ;
00875             qmLSSResSize = molPtr->get_qmLSSResSize() ;
00876             qmLSSRefIDs = molPtr->get_qmLSSRefIDs() ;
00877             qmLSSRefSize = molPtr->get_qmLSSRefSize() ;
00878             
00879             lssPrepare();
00880         }
00881         
00882         numArrivedQMMsg = 0 ;
00883         numArrivedPntChrgMsg = 0 ;
00884         
00885         qmCoord = new ComputeQMAtom[numQMAtms];
00886         
00887         replaceForces = 0;
00888         if (molPtr->get_qmReplaceAll()) {
00889             replaceForces = 1;
00890         }
00891         
00892         cSMDon = molPtr->get_qmcSMD() ;
00893         if (cSMDon) {
00894                 
00895                 // We have to initialize the guide particles during the first step.
00896                 cSMDInitGuides = 0;
00897                 
00898                 cSMDnumInstances = molPtr->get_cSMDnumInst();
00899                 cSMDindex = molPtr->get_cSMDindex();
00900                 cSMDindxLen = molPtr->get_cSMDindxLen();
00901                 cSMDpairs = molPtr->get_cSMDpairs(); 
00902                 cSMDKs = molPtr->get_cSMDKs();
00903                 cSMDVels = molPtr->get_cSMDVels();
00904                 cSMDcoffs = molPtr->get_cSMDcoffs();
00905                 
00906                 cSMDguides = new Position[cSMDnumInstances];
00907                 cSMDForces = new Force[cSMDnumInstances];
00908         }
00909                 
00910                 
00911         DebugM(4,"Initializing DCD file for charge information." << std::endl);
00912         
00913         // Initializes output DCD file for charge information.
00914         if (simParams->qmOutFreq > 0) {
00915             std::string dcdOutFileStr;
00916             dcdOutFileStr.clear();
00917             dcdOutFileStr.append(simParams->outputFilename) ;
00918             dcdOutFileStr.append(".qdcd") ;
00919             dcdOutFile = open_dcd_write(dcdOutFileStr.c_str()) ;
00920             
00921             if (dcdOutFile == DCD_FILEEXISTS) {
00922                 iout << iERROR << "DCD file " << dcdOutFile << " already exists!!\n" << endi;
00923                 NAMD_err("Could not write QM charge DCD file.");
00924             } else if (dcdOutFile < 0) {
00925                 iout << iERROR << "Couldn't open DCD file " << dcdOutFile << ".\n" << endi;
00926                 NAMD_err("Could not write QM charge DCD file.");
00927             } else if (! dcdOutFile) {
00928                 NAMD_bug("ComputeQMMgr::recvPartQM open_dcd_write returned fileid of zero");
00929             }
00930             
00931             timeStep = simParams->firstTimestep;
00932             
00933             int NSAVC, NFILE, NPRIV, NSTEP;
00934             NSAVC = simParams->qmOutFreq;
00935             NPRIV = timeStep;
00936             NSTEP = NPRIV - NSAVC;
00937             NFILE = 0;
00938             
00939             //  Write out the header
00940             int ret_code = write_dcdheader(dcdOutFile, dcdOutFileStr.c_str(),
00941             numQMAtms, NFILE, NPRIV, NSAVC, NSTEP,
00942             simParams->dt/TIMEFACTOR, 0);
00943 
00944             if (ret_code<0) {
00945                 NAMD_err("Writing of DCD header failed!!");
00946             }
00947             
00948             // The DCD write function takes 3 independent arrays for X, Y and Z
00949             // coordinates, but we allocate one and send in the pieces.
00950             outputData = new Real[3*numQMAtms];
00951         }
00952         
00953         DebugM(4,"Initializing DCD file for position information." << std::endl);
00954         // Initializes output DCD file for position information.
00955         if (simParams->qmPosOutFreq > 0) {
00956             std::string dcdPosOutFileStr;
00957             dcdPosOutFileStr.clear();
00958             dcdPosOutFileStr.append(simParams->outputFilename) ;
00959             dcdPosOutFileStr.append(".QMonly.dcd") ;
00960             dcdPosOutFile = open_dcd_write(dcdPosOutFileStr.c_str()) ;
00961             
00962             if (dcdPosOutFile == DCD_FILEEXISTS) {
00963                 iout << iERROR << "DCD file " << dcdPosOutFile << " already exists!!\n" << endi;
00964                 NAMD_err("Could not write QM charge DCD file.");
00965             } else if (dcdPosOutFile < 0) {
00966                 iout << iERROR << "Couldn't open DCD file " << dcdPosOutFile << ".\n" << endi;
00967                 NAMD_err("Could not write QM charge DCD file.");
00968             } else if (! dcdPosOutFile) {
00969                 NAMD_bug("ComputeQMMgr::recvPartQM open_dcd_write returned fileid of zero");
00970             }
00971             
00972             timeStep = simParams->firstTimestep;
00973             
00974             int NSAVC, NFILE, NPRIV, NSTEP;
00975             NSAVC = simParams->qmOutFreq;
00976             NPRIV = timeStep;
00977             NSTEP = NPRIV - NSAVC;
00978             NFILE = 0;
00979             
00980             //  Write out the header
00981             int ret_code = write_dcdheader(dcdPosOutFile, dcdPosOutFileStr.c_str(),
00982             numQMAtms, NFILE, NPRIV, NSAVC, NSTEP,
00983             simParams->dt/TIMEFACTOR, 0);
00984 
00985             if (ret_code<0) {
00986                 NAMD_err("Writing of DCD header failed!!");
00987             }
00988             
00989             // The DCD write function takes 3 independent arrays for X, Y and Z
00990             // coordinates, but we allocate one and send in the pieces.
00991             outputData = new Real[3*numQMAtms];
00992         }
00993         
00994         // Prepares list of PEs which will run the QM software
00995         int simsPerNode = simParams->qmSimsPerNode ;
00996         int zeroNodeSize = CmiNumPesOnPhysicalNode(0);
00997         
00998         // Check if the node has enought PEs to run the requested number of simulations.
00999         if ( zeroNodeSize < simsPerNode ) {
01000             iout << iERROR << "There are " << zeroNodeSize << " cores pernode, but "
01001             << simsPerNode << " QM simulations per node were requested.\n" << endi ;
01002             NAMD_die("Error preparing QM simulations.");
01003         }
01004         
01005         DebugM(4,"Preparing PE list for QM execution.\n");
01006         qmPEs.clear(); // Making sure its empty.
01007         
01008         int numNodes = CmiNumPhysicalNodes();
01009         int numPlacedQMGrps = 0;
01010         int placedOnNode = 0;
01011         int nodeIt = 0 ;
01012         
01013         // The default is to only run on rank zero.
01014         if ( simsPerNode <= 0 ) {
01015             qmPEs.push_back(0);
01016             numPlacedQMGrps = 1;
01017         }
01018         
01019         while ( (numPlacedQMGrps < numQMGrps) && (simsPerNode > 0) ) {
01020             
01021             // If we searched all nodes, break the loop.
01022             if (nodeIt == numNodes) {
01023                 break;
01024             }
01025             
01026             int *pelist;
01027             int nodeSize;
01028             CmiGetPesOnPhysicalNode(CmiPhysicalNodeID(nodeIt), &pelist, &nodeSize);
01029             
01030             DebugM(4,"Checking node " << nodeIt +1 << " out of " << numNodes
01031             << " (" << nodeSize << " PEs: " << pelist[0] << " to " 
01032             << pelist[nodeSize-1] << ")." << std::endl );
01033             
01034             for ( int i=0; i<nodeSize; ++i ) {
01035                 
01036                 // Check if the PE has patches. We only run on PEs with patches!
01037                 if ( (PatchMap::Object())->numPatchesOnNode(pelist[i]) == 0 ) {
01038                     DebugM(1,"PE " << pelist[i] << " has no patches! We'll skip it." 
01039                     << std::endl);
01040                     continue ;
01041                 }
01042                 
01043                 // Add the target PE on the target node to the list
01044                 // of PEs which will carry QM simulations.
01045                 qmPEs.push_back(pelist[i]);
01046                 
01047                 DebugM(1,"Added PE to QM execution list: " << pelist[i]  << "\n");
01048                 
01049                 numPlacedQMGrps++;
01050                 placedOnNode++;
01051                 
01052                 if (placedOnNode == simsPerNode) {
01053                     DebugM(1,"Placed enought simulations on this node.\n");
01054                     break;
01055                 }
01056                 
01057                 
01058             }
01059             
01060             nodeIt++ ;
01061             placedOnNode = 0;
01062         }
01063         
01064         if ( numPlacedQMGrps < numQMGrps ) {
01065             iout << iWARN << "Could not compute all QM groups in parallel.\n" << endi ;
01066         }
01067         
01068         iout << iINFO << "List of ranks running QM simulations: " << qmPEs[0] ;
01069         for (int i=1; i < qmPEs.size(); i++) {
01070             iout << ", " << qmPEs[i] ;
01071         }
01072         iout << ".\n" << endi;
01073         
01074     }
01075     
01076     DebugM(1,"Receiving from rank " << msg->sourceNode
01077     << " a total of " << msg->numAtoms << " QM atoms." << std::endl);
01078     
01079     // In case we are NOT using point charges but there are QM-MM bonds,
01080     // test each QM message for MM1 atoms.
01081     if (meNumMMIndx > 0) {
01082         
01083         ResizeArray< ComputeQMAtom > tmpAtm;
01084         ComputeQMAtom *data_ptr = msg->coord;
01085         
01086         for (int i=0; i<msg->numAtoms; i++) {
01087             if (data_ptr[i].vdwType < 0) {
01088                 tmpAtm.add(data_ptr[i]) ;
01089             }
01090         }
01091         
01092         QMPntChrgMsg *pntChrgMsg = new (tmpAtm.size(), 0) QMPntChrgMsg;
01093         pntChrgMsg->sourceNode = msg->sourceNode ;
01094         pntChrgMsg->numAtoms = tmpAtm.size() ;
01095         ComputeQMPntChrg* newPCData = pntChrgMsg->coord ;
01096         
01097         QMCoordMsg *newMsg = msg;
01098         
01099         if (tmpAtm.size() > 0) {
01100             
01101             newMsg = new (msg->numAtoms - tmpAtm.size(),0, 0) QMCoordMsg;
01102             newMsg->sourceNode = msg->sourceNode ;
01103             newMsg->numAtoms = msg->numAtoms - tmpAtm.size() ;
01104             newMsg->timestep = msg->timestep ;
01105             ComputeQMAtom *newMsgData = newMsg->coord;
01106             
01107             for (int i=0; i<msg->numAtoms; i++) {
01108                 if (data_ptr[i].vdwType < 0) {
01109                     newPCData->position = data_ptr[i].position ;
01110                     newPCData->charge = data_ptr[i].charge ;
01111                     newPCData->id = data_ptr[i].id ;
01112                     newPCData->qmGrpID = data_ptr[i].qmGrpID ;
01113                     newPCData->homeIndx = data_ptr[i].homeIndx ;
01114                     newPCData->dist = 0 ;
01115                     newPCData->mass = 0 ;
01116                     newPCData->vdwType = 0 ;
01117                     newPCData++;
01118                 }
01119                 else {
01120                     *newMsgData = data_ptr[i] ;
01121                     newMsgData++;
01122                 }
01123             }
01124             
01125             delete msg;
01126             
01127         }
01128         
01129         qmCoordMsgs[numArrivedQMMsg] = newMsg;
01130         ++numArrivedQMMsg;
01131         
01132         pntChrgCoordMsgs[numArrivedPntChrgMsg] = pntChrgMsg;
01133         ++numArrivedPntChrgMsg;
01134     }
01135     else {
01136         qmCoordMsgs[numArrivedQMMsg] = msg;
01137         ++numArrivedQMMsg;
01138     }
01139     
01140     if ( numArrivedQMMsg < numSources ) 
01141         return;
01142     
01143     // Now that all messages arrived, get all QM positions.
01144     for (int msgIt=0; msgIt < numArrivedQMMsg; msgIt++){
01145         
01146         DebugM(1, "Getting positions for " << qmCoordMsgs[msgIt]->numAtoms 
01147         << " QM atoms in this message." << std::endl);
01148         
01149         for ( int i=0; i < qmCoordMsgs[msgIt]->numAtoms; ++i ) {
01150             qmCoord[qmAtmIter] = qmCoordMsgs[msgIt]->coord[i];
01151             qmAtmIter++;
01152         }
01153     }
01154     
01155     if (qmAtmIter != numQMAtms) {
01156         iout << iERROR << "The number of QM atoms received (" << qmAtmIter
01157         << ") is different than expected: " << numQMAtms << "\n" << endi;
01158         NAMD_err("Problems broadcasting QM atoms.");
01159     }
01160     
01161     // Resets the counter for the next step.
01162     numArrivedQMMsg = 0;
01163     qmAtmIter = 0;
01164     
01165     timeStep = qmCoordMsgs[0]->timestep;
01166     
01167     // Makes sure there is no trash or old info in the force array.
01168     // This is where we accumulate forces from the QM software and our own
01169     // Coulomb forces. It will have info on QM atoms and point charges only.
01170     for (int i=0; i<numAtoms; ++i ) {
01171         force[i].force = 0;  // this assigns 0 (ZERO) to all componenets of the vector.
01172         force[i].replace = 0;  // By default, no atom has all forces replaced.
01173         force[i].homeIndx = -1;  // To prevent errors from sliping.
01174         force[i].charge = 0;
01175         force[i].id = i;        // Initializes the variable since not all forces are sent to all ranks.
01176     }
01177     
01178     
01179     for (int i=0; i<numQMAtms; i++) {
01180         // Each force receives the home index of its atom with respect to the 
01181         // local set of atoms in each node.
01182         if (force[qmCoord[i].id].homeIndx != -1
01183             && force[qmCoord[i].id].homeIndx != qmCoord[i].homeIndx
01184         ) {
01185             iout << iERROR << "Overloading QM atom " 
01186             << qmCoord[i].id << "; home index: " 
01187             << force[qmCoord[i].id].homeIndx << " with " << qmCoord[i].homeIndx
01188             << "\n" << endi ;
01189             NAMD_die("Error preparing QM simulations.");
01190         }
01191         
01192         force[qmCoord[i].id].homeIndx = qmCoord[i].homeIndx;
01193         // Each force on QM atoms has an indicator so NAMD knows if all 
01194         // NAMD forces should be erased and completely overwritten by the
01195         // external QM forces.
01196         force[qmCoord[i].id].replace = replaceForces;
01197     }
01198     
01199     if (noPC) {
01200         // this pointer should be freed in rank zero, after receiving it.
01201         QMPntChrgMsg *pntChrgMsg = new (0, 0) QMPntChrgMsg;
01202         pntChrgMsg->sourceNode = CkMyPe();
01203         pntChrgMsg->numAtoms = 0;
01204         
01205         QMProxy[0].recvPntChrg(pntChrgMsg);
01206     }
01207     else {
01208         // The default mode is to update the poitn charge selection at every step.
01209         int pcSelMode = PCMODEUPDATESEL;
01210         
01211         int msgPCListSize = 0;
01212         // We check wether point charges are to be selected with a stride.
01213         if ( qmPCFreq > 0 ) {
01214             if (timeStep % qmPCFreq == 0) {
01215                 // Marks that the PC list determined in this step will
01216                 // be broadcast on the *next* step, and kept for the following
01217                 // qmPCFreq-1 steps.
01218                 resendPCList = true;
01219                 
01220                 // Clears the list since we don't know how many charges 
01221                 // will be selected this time.
01222                 delete [] pcIDList;
01223             }
01224             else {
01225                 // If the PC selection is not to be updated in this step,
01226                 // inform the nodes that the previous list (or the updated
01227                 // list being sent in this message) is to be used and only 
01228                 // updated positions will be returned.
01229                 pcSelMode = PCMODEUPDATEPOS;
01230                 
01231                 // If this is the first step after a PC re-selection, all 
01232                 // ranks receive the total list, since atoms may move between
01233                 // PEs in between PC re-assignments (if they are far enought apart
01234                 // or if you are unlucky)
01235                 if (resendPCList) {
01236                     msgPCListSize = pcIDListSize;
01237                     resendPCList = false;
01238                 }
01239             }
01240         }
01241         
01242         // In case we are using custom selection of point charges, indicate the mode.
01243         if (simParams->qmCustomPCSel)
01244             pcSelMode = PCMODECUSTOMSEL;
01245         
01246         DebugM(1,"Broadcasting current positions of QM atoms.\n");
01247         for ( int j=0; j < numSources; ++j ) {
01248             // This pointer will be freed in the receiving rank.
01249             QMCoordMsg *qmFullMsg = new (numQMAtms, msgPCListSize, 0) QMCoordMsg;
01250             qmFullMsg->sourceNode = CkMyPe();
01251             qmFullMsg->numAtoms = numQMAtms;
01252             qmFullMsg->pcSelMode = pcSelMode;
01253             qmFullMsg->numPCIndxs =  msgPCListSize;
01254             ComputeQMAtom *data_ptr = qmFullMsg->coord;
01255             int *msgPCListP = qmFullMsg->pcIndxs;
01256             
01257             for (int i=0; i<numQMAtms; i++) {
01258                 data_ptr->position = qmCoord[i].position;
01259                 data_ptr->charge = qmCoord[i].charge;
01260                 data_ptr->id = qmCoord[i].id;
01261                 data_ptr->qmGrpID = qmCoord[i].qmGrpID;
01262                 data_ptr->homeIndx = -1; // So there is no mistake that there is no info here.
01263                 data_ptr++;
01264             }
01265             
01266             for (int i=0; i<msgPCListSize; i++) {
01267                 msgPCListP[i] = pcIDList[i] ;
01268             }
01269             
01270             #ifdef DEBUG_QM
01271             if (j == 0)
01272                 Write_PDB("qmMsg.pdb", qmFullMsg) ;
01273             #endif
01274             
01275             // The messages are deleted later, we will need them.
01276             QMProxy[qmCoordMsgs[j]->sourceNode].recvFullQM(qmFullMsg);
01277         }
01278     }
01279 }
01280 
01281 void ComputeQMMgr::recvFullQM(QMCoordMsg* qmFullMsg) {
01282     
01283     if (subsArray.size() > 0)
01284         subsArray.clear();
01285     
01286     QMCompute->processFullQM(qmFullMsg);
01287 }
01288 
01289 typedef std::map<Real, std::pair<Position,BigReal> > GrpDistMap ;
01290 typedef std::pair<Position,BigReal> PosDistPair ;
01291 
01292 void ComputeQM::processFullQM(QMCoordMsg* qmFullMsg) {
01293     
01294     ResizeArrayIter<PatchElem> ap(patchList); 
01295     
01296     const Lattice & lattice = patchList[0].p->lattice;
01297     
01298     // Dynamically accumulates point charges as they are found.
01299     // The same MM atom may be added to this vector more than once if it is 
01300     // within the cutoff region of two or more QM regions.
01301     ResizeArray<ComputeQMPntChrg> compPCVec ;
01302     
01303     // This will keep the set of QM groups with which each 
01304     // point charge should be used. It is re-used for each atom in this
01305     // patch so we can controll the creation of the compPCVec vector.
01306     // The first item in the pair is the QM Group ID, the second is a pair 
01307     // with the point charge position (after PBC wraping) and the shortest
01308     // distance between the point charge and an atom of this QM group.
01309     GrpDistMap chrgGrpMap ;
01310     
01311     DebugM(4,"Rank " << CkMyPe() << " receiving from rank " << qmFullMsg->sourceNode
01312     << " a total of " << qmFullMsg->numAtoms << " QM atoms and " 
01313     << qmFullMsg->numPCIndxs << " PC IDs." << std::endl);
01314     
01315     // If this is the firts step after a step of PC re-selection,
01316     // we store all PC IDs that all PEs found, in case some atoms end up
01317     // here untill the next re-selection step.
01318     if (qmFullMsg->numPCIndxs) {
01319         
01320         pcIDSortList.clear();
01321         
01322         int *pcIndx = qmFullMsg->pcIndxs;
01323         for (int i=0; i< qmFullMsg->numPCIndxs;i++) {
01324             pcIDSortList.load(pcIndx[i]);
01325         }
01326         
01327         pcIDSortList.sort();
01328     }
01329     
01330     int totalNodeAtoms = 0;
01331     int atomIter = 0;
01332     int uniqueCounter = 0;
01333     
01334     const ComputeQMAtom *qmDataPtn = qmFullMsg->coord;
01335     
01336     switch ( qmFullMsg->pcSelMode ) {
01337     
01338     case PCMODEUPDATESEL:
01339     {
01340         
01341     DebugM(4,"Updating PC selection.\n")
01342     
01343     // Loops over all atoms in this node and checks if any MM atom is within 
01344     // the cutof radius from a QM atom
01345     for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
01346         
01347         CompAtom *x = (*pdIt).compAtomP;
01348         const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
01349         int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
01350         const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
01351         
01352         totalNodeAtoms += localNumAtoms;
01353         
01354         // Iterates over the local atoms in a PE, all patches.
01355         for(int i=0; i<localNumAtoms; ++i) {
01356             
01357             // A new subset for each atom in this node.
01358             chrgGrpMap.clear();
01359             
01360             Real pcGrpID = qmAtomGroup[xExt[i].id];
01361             
01362             Real charge = x[i].charge;
01363             
01364             // If the atom is a QM atom, there will be no charge info in the 
01365             // atom box. Therefore, we take the charge in the previous step
01366             // in case this atom is a point charge for a neighboring QM region.
01367             if (pcGrpID > 0) {
01368                 for (int qi=0; qi<numQMAtms; qi++) {
01369                     if (qmAtmIndx[qi] == xExt[i].id) {
01370                         charge = qmAtmChrg[qi];
01371                         break;
01372                     }
01373                 }
01374             }
01375             
01376             qmDataPtn = qmFullMsg->coord;
01377             for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
01378                 
01379                 Real qmGrpID = qmDataPtn->qmGrpID;
01380                 
01381                 // We check If a QM atom and the node atom "i"
01382                 // belong to the same QM group. If not, atom "i" is either an MM
01383                 // atom or a QM atom from another group, which will be seen
01384                 // as an MM point charge.
01385                 // If they belong to the same group, just skip the distance 
01386                 // calculation and move on to the next QM atom.
01387                 // This loop needs to go over all QM atoms since a point charge
01388                 // may me be close to two different QM groups, in which case it 
01389                 // will be sent to both.
01390                 if (qmGrpID == pcGrpID) {
01391                     continue;
01392                 }
01393                 
01394                 // calculate distance between QM atom and Poitn Charge
01395                 Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
01396                 
01397                 BigReal dist = r12.length();  // Distance between atoms
01398                 
01399                 if ( dist < cutoff ){
01400                     
01401                     // Since the QM region may have been wrapped around the periodic box, we
01402                     // reposition point charges with respect to unwrapped coordinates of QM atoms,
01403                     // finding the shortest distance between the point charge and the QM region.
01404                     
01405                     Position posMM = qmDataPtn->position + r12;
01406                     auto ret = chrgGrpMap.find(qmGrpID) ;
01407                     
01408                     // If 'ret' points to end, it means the item was not found
01409                     // and will be added, which means a new QM region has this 
01410                     // atom within its cutoff region,
01411                     // which means a new point charge will be 
01412                     // created in the QM system in question.
01413                     // 'ret' means a lot to me!
01414                     if ( ret ==  chrgGrpMap.end()) {
01415                         chrgGrpMap.insert(std::pair<Real,PosDistPair>(qmGrpID, 
01416                                                          PosDistPair(posMM,dist) ) );
01417                     }
01418                     else {
01419                         // If the current distance is smaller than the one in 
01420                         // the map, we update it so that we have the smallest
01421                         // distance between the point charge and any QM atom in 
01422                         // this group.
01423                         if (dist < ret->second.second) {
01424                             ret->second.first = posMM ;
01425                             ret->second.second = dist ;
01426                         }
01427                     }
01428                 }
01429             }
01430             
01431             for(auto mapIt = chrgGrpMap.begin(); 
01432                 mapIt != chrgGrpMap.end(); mapIt++) {
01433                 
01434                 // We now add the final info about this point charge 
01435                 // to the vector, repeating it for each QM group that has it
01436                 // within its cuttoff radius.
01437                 compPCVec.add(
01438                     ComputeQMPntChrg(mapIt->second.first, charge, xExt[i].id,
01439                                   mapIt->first, atomIter, mapIt->second.second, 
01440                                   fullAtms[i].mass, fullAtms[i].vdwType)
01441                                );
01442             }
01443             
01444             // Counts how many atoms are seens as point charges, by one or more
01445             // QM groups.
01446             if (chrgGrpMap.size() > 0)
01447                 uniqueCounter++;
01448             
01449             atomIter++;
01450         }
01451         
01452         (*pdIt).posBoxP->close(&x);
01453     }
01454     
01455     }break; // End case PCMODEUPDATESEL
01456     
01457     case PCMODEUPDATEPOS: 
01458     {
01459     
01460     DebugM(4,"Updating PC positions ONLY.\n")
01461     
01462     for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
01463         
01464         CompAtom *x = (*pdIt).compAtomP;
01465         const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
01466         int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
01467         const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
01468         
01469         totalNodeAtoms += localNumAtoms;
01470         
01471         // Iterates over the local atoms in a PE, all patches.
01472         for(int i=0; i<localNumAtoms; ++i) {
01473             
01474             if (pcIDSortList.find(xExt[i].id) != NULL ) {
01475                 
01476                 BigReal dist =  0;
01477                 Position posMM = 0; // Temporary point charge position for QM calculation.
01478                 
01479                 bool firstIngroupQM = true;
01480                 qmDataPtn = qmFullMsg->coord;
01481                 for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
01482                     
01483                     // Assuming we only have one QM region for now
01484                     // This optio will be deprecated
01485                     
01486                     // Only verify distances to atoms in the QM group for which we are
01487                     // gathering point charges.
01488 //                     if ( qmDataPtn->qmGrpID != qmGrpIDArray[grpIndx] ) continue;
01489                     
01490                     // calculate distance between QM atom and Poitn Charge
01491                     Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
01492                     
01493                     // We wait to find the first QM atom in the target QM group before we
01494                     // set the first guess at a minimal distance and point charge position.
01495                     if ( firstIngroupQM ) {
01496                         dist = r12.length();
01497                         
01498                         // Reposition classical point charge with respect to unwrapped QM position.
01499                         posMM = qmDataPtn->position + r12;
01500                         
01501                         firstIngroupQM = false;
01502                     }
01503                     
01504                     // If we find a QM atom with a shorter dstance to the point charge,
01505                     // set that as the PC distance and reposition the PC
01506                     if ( r12.length() < dist ) {
01507                         dist = r12.length();
01508                         posMM = qmDataPtn->position + r12 ;
01509                     }
01510                     
01511                 }
01512                 
01513                 Real pcGrpID = qmAtomGroup[xExt[i].id];
01514                 Real charge = x[i].charge;
01515                 
01516                 // If the atom is a QM atom, there will be no charge info in the 
01517                 // atom box. Therefore, we take the charge in the previous step
01518                 // in case this atom is a point charge for a neighboring QM region.
01519                 if (pcGrpID > 0) {
01520                     for (int qi=0; qi<numQMAtms; qi++) {
01521                         if (qmAtmIndx[qi] == xExt[i].id) {
01522                             charge = qmAtmChrg[qi];
01523                             break;
01524                         }
01525                         
01526                     }
01527                 }
01528                 
01529                 compPCVec.add(
01530                     ComputeQMPntChrg(posMM, charge, xExt[i].id,
01531                                   0, atomIter, 0, 
01532                                   fullAtms[i].mass, fullAtms[i].vdwType));
01533             }
01534             
01535             atomIter++;
01536         }
01537         
01538         (*pdIt).posBoxP->close(&x);
01539     }
01540     }break ; // End case PCMODEUPDATEPOS
01541     
01542     case PCMODECUSTOMSEL:
01543     {
01544     
01545     DebugM(4,"Updating PC positions for custom PC selection.\n")
01546     
01547     for (auto pdIt = patchData.begin(); pdIt != patchData.end(); pdIt++ ) {
01548         
01549         CompAtom *x = (*pdIt).compAtomP;
01550         const CompAtomExt *xExt =  (*pdIt).homePatchP->getCompAtomExtInfo();
01551         int localNumAtoms =  (*pdIt).homePatchP->getNumAtoms();
01552         const FullAtom* fullAtms = (*pdIt).homePatchP->getAtomList().const_begin() ;
01553         
01554         totalNodeAtoms += localNumAtoms;
01555         
01556         // Iterates over the local atoms in a PE, all patches.
01557         for(int i=0; i<localNumAtoms; ++i) {
01558             
01559             // Checks if the atom ID belongs to a point charge list,
01560             // which indicates it is to be sent to rank zero for QM calculations.
01561             // With one list per QM group, the same point charge can be added to 
01562             // different QM groups, as if it was being dynamically selected.
01563             for (int grpIndx=0; grpIndx<numQMGrps; grpIndx++) {
01564                 
01565                 if (customPCLists[grpIndx].find(xExt[i].id) != NULL){
01566                     
01567                     // Initialize the search for the smallest distance between
01568                     // point charge and QM atoms. This will determine the closest position
01569                     // of the point charge accounting for periodic boundary conditions.
01570                     
01571                     // Since the QM region may have been wrapped around the periodic box, we
01572                     // reposition point charges with respect to unwrapped coordinates of QM atoms.
01573                     
01574                     BigReal dist =  0;
01575                     Position posMM = 0; // Temporary point charge position for QM calculation.
01576                     
01577                     bool firstIngroupQM = true;
01578                     qmDataPtn = qmFullMsg->coord;
01579                     for(int j=0; j<qmFullMsg->numAtoms; j++, ++qmDataPtn) {
01580                         
01581                         // Only verify distances to atoms in the QM group for which we are
01582                         // gathering point charges.
01583                         if ( qmDataPtn->qmGrpID != qmGrpIDArray[grpIndx] ) continue;
01584                         
01585                         // calculate distance between QM atom and Poitn Charge
01586                         Vector r12 = lattice.delta(x[i].position,qmDataPtn->position);
01587                         
01588                         // We wait to find the first QM atom in the target QM group before we
01589                         // set the first guess at a minimal distance and point charge position.
01590                         if ( firstIngroupQM ) {
01591                             dist = r12.length();
01592                             
01593                             // Reposition classical point charge with respect to unwrapped QM position.
01594                             posMM = qmDataPtn->position + r12;
01595                             
01596                             firstIngroupQM = false;
01597                         }
01598                         
01599                         // If we find a QM atom with a shorter dstance to the point charge,
01600                         // set that as the PC distance and reposition the PC
01601                         if ( r12.length() < dist ) {
01602                             dist = r12.length();
01603                             posMM = qmDataPtn->position + r12 ;
01604                         }
01605                         
01606                     }
01607                     
01608                     // After determining the position o fthe point charge, get its partial charge.
01609                     Real pcGrpID = qmAtomGroup[xExt[i].id];
01610                     Real charge = x[i].charge;
01611                     
01612                     // If the atom is a QM atom, there will be no charge info in the 
01613                     // atom box. Therefore, we take the charge in the previous step
01614                     // in case this atom is a point charge for a neighboring QM region.
01615                     if (pcGrpID > 0) {
01616                         for (int qi=0; qi<numQMAtms; qi++) {
01617                             if (qmAtmIndx[qi] == xExt[i].id) {
01618                                 charge = qmAtmChrg[qi];
01619                                 break;
01620                             }
01621                         }
01622                     }
01623                     
01624                     compPCVec.add(
01625                         ComputeQMPntChrg(posMM, charge, xExt[i].id,
01626                                       qmGrpIDArray[grpIndx], atomIter, dist, 
01627                                       fullAtms[i].mass, fullAtms[i].vdwType));
01628                 }
01629                 
01630             }
01631             
01632             atomIter++;
01633         }
01634         
01635         (*pdIt).posBoxP->close(&x);
01636     }
01637     
01638     } break;
01639     }
01640     
01641     DebugM(4,"Rank " << CkMyPe() << " found a total of " << compPCVec.size()
01642     << " point charges, out of " << totalNodeAtoms
01643     << " atoms in this node. " << uniqueCounter << " are unique." << std::endl);
01644     
01645     // Send only the MM atoms within radius
01646     
01647     // this pointer should be freed in rank zero, after receiving it.
01648     QMPntChrgMsg *pntChrgMsg = new (compPCVec.size(), 0) QMPntChrgMsg;
01649     pntChrgMsg->sourceNode = CkMyPe();
01650     pntChrgMsg->numAtoms = compPCVec.size();
01651     
01652     for (int i=0; i<compPCVec.size(); i++ ) {
01653         
01654         // Only sends the positions and charges of atoms within
01655         // cutoff of a (any) QM atom. Repeats for the number of QM groups
01656         // this charge is near to, and indicates the QM group it should 
01657         // be used with.
01658         pntChrgMsg->coord[i] = compPCVec[i] ;
01659         
01660     }
01661     
01662     DebugM(4,"Rank " << pntChrgMsg->sourceNode << " sending a total of " 
01663     << compPCVec.size() << " elements to rank zero." << std::endl);
01664     
01665     CProxy_ComputeQMMgr QMProxy(CkpvAccess(BOCclass_group).computeQMMgr);
01666     QMProxy[0].recvPntChrg(pntChrgMsg);
01667     
01668     delete qmFullMsg;
01669 }
01670 
01671 
01672 
01673 void ComputeQMMgr::procBonds(int numBonds,
01674                              const int *const qmGrpBondsGrp,
01675                              const int *const qmMMBondedIndxGrp,
01676                              const int *const *const chargeTarget,
01677                              const int *const numTargs,
01678                              const QMPCVec grpPntChrgVec,
01679                              QMPCVec &grpAppldChrgVec) {
01680     
01681     DebugM(1,"Processing QM-MM bonds in rank zero.\n");
01682     
01683     // Indices and IDs of MM atoms which will be passed as point charges.
01684     std::map<int, int> mmPntChrgMap ;
01685     
01686     // Build the map of atom ID and charge for this QM group's point charges.
01687     // They may be changed by different charge distribution schemes and/or
01688     // dummy aotm placements.
01689     for (size_t i=0; i<grpPntChrgVec.size(); i++) {
01690         
01691         mmPntChrgMap.insert(std::pair<int,int>(grpPntChrgVec[i].id, (int) i) );
01692         
01693         grpAppldChrgVec.push_back( grpPntChrgVec[i] ) ;
01694         
01695     }
01696     
01697     // If we are treating QM-MM bonds and if there are any
01698     // in this QM group, we have to treat the MM partial charge
01699     // using some scheme.
01700     for (int bondIt = 0; bondIt < numBonds; bondIt++) {
01701         
01702         // Gets the global index of this particular QM-MM bond.
01703         int bondIndx = qmGrpBondsGrp[bondIt] ;
01704         
01705         auto retIt = mmPntChrgMap.find(qmMMBondedIndxGrp[bondIt]) ;
01706         
01707         // Checks if this MM atom was included as a 
01708         // point charge due proximity to a QM atom.
01709         if (retIt == mmPntChrgMap.end()) {
01710             // If it wasn't, there is an error somwhere.
01711             
01712             iout << iERROR << "The MM atom " << qmMMBondedIndxGrp[bondIt]
01713             << " is bound to a QM atom, but it was not selected as a poitn charge."
01714             << " Check your cutoff radius!\n" << endi ;
01715             
01716             NAMD_die("Charge placement error in QM-MM bond.");
01717         }
01718         
01719         // Gets the (local) index of the MM atom in the point charge vector.
01720         int mmIndex = (*retIt).second;
01721         // Gets the position of the MM atom and its partial charge.
01722         Position mmPos = grpAppldChrgVec[mmIndex].position ;
01723         BigReal mmCharge = grpAppldChrgVec[mmIndex].charge/numTargs[bondIndx] ;
01724         
01725         // gives part of the MM charge to neighboring atoms
01726         for (int i=0; i<numTargs[bondIndx]; i++){
01727             
01728             int targetIndxGLobal = chargeTarget[bondIndx][i] ;
01729             
01730             retIt = mmPntChrgMap.find(targetIndxGLobal);
01731             
01732             // If one of the neighboring atoms which should receive charge
01733             // is not among partial charges, stop and run away.
01734             if (retIt == mmPntChrgMap.end()) {
01735                 
01736                 iout << iERROR << "The MM atom " << targetIndxGLobal
01737                 << " is bound to the MM atom of a QM-MM bond and is needed for"
01738                 << " the required bond scheme"
01739                 << " but it was not selected as a poitn charge."
01740                 << " Check your cutoff radius!\n" << endi ;
01741                 
01742                 NAMD_die("Charge placement error in QM-MM bond.");
01743             }
01744             
01745             int trgIndxLocal = (*retIt).second;
01746             
01747             // Choose charge treatment method and apply it.
01748             switch (simParams->qmBondScheme) {
01749             
01750                 // Charge Schiftig scheme.
01751                 case QMSCHEMECS:
01752                 {
01753                 
01754                 // Charge Shifting Scheme (DOI: 10.1021/ct100530r)
01755                 // Here we create a dipole to counter act the charge movement
01756                 // we created by moving parts of the MM charge to target MM atoms.
01757                 // 
01758                 // The MM1 charge is equally divided and added to all MM2 atoms.
01759                 // Two point charges are created in the direction of the MM1-MM2 bond,
01760                 // one before and one after the MM2 atom.
01761                 // 
01762                 // Below we see the diagram of the positions of new charges along
01763                 // the direction of the bond between the MM1 atom and a target 
01764                 // MM2 atom, wiht respect to the new Dummy atom (a Hydrogen).
01765                 // 
01766                 //  QM --------  MM1(p0) ------------ MM2
01767                 //  QM ------ H ------------ (+p0) - (MM2 +p0) - (-p0)
01768                 
01769                 Position trgPos = grpAppldChrgVec[trgIndxLocal].position ;
01770                 
01771                 // We now store in MM2 to which PC it is bound (the index for MM1).
01772                 // This will be used to decompose virtual PC forces onto MM1 and MM2.
01773                 grpAppldChrgVec[trgIndxLocal].qmGrpID = -2;
01774                 grpAppldChrgVec[trgIndxLocal].homeIndx = mmIndex;
01775                 
01776                 // We create a new point charge at the same position so that
01777                 // the fraction of charge from the MM1 atom is "added" to the
01778                 // MM2 position, but without actually changing the charge of 
01779                 // the original MM2 atom. This trick helps keeping the original 
01780                 // charge for electrostatic calculations after the new charges 
01781                 // of QM atoms is received from the QM software.
01782                 Real Cq0 = 1.0;
01783                 grpAppldChrgVec.push_back(
01784                     // Cq is stored in the distance slot, since it has not been computed
01785                     // for a virtual charge.
01786                     ComputeQMPntChrg(trgPos, mmCharge, trgIndxLocal, 0, -1, Cq0, 0, 0)
01787                 );
01788                 
01789                 Vector bondVec = trgPos - mmPos ;
01790                 
01791                 // For Cq-plus virtual charge
01792                 Real Cqp = 0.94;
01793                 // For Cq-minus virtual charge
01794                 Real Cqm = 1.06;
01795                 Vector bondVec1 = bondVec*Cqp ;
01796                 Vector bondVec2 = bondVec*Cqm ;
01797                 
01798                 Position chrgPos1 = mmPos + bondVec1;
01799                 Position chrgPos2 = mmPos + bondVec2;
01800                 
01801                 BigReal trgChrg1 = mmCharge;
01802                 BigReal trgChrg2 = -1*mmCharge;
01803                 
01804                 grpAppldChrgVec.push_back(
01805                     ComputeQMPntChrg(chrgPos1, trgChrg1, trgIndxLocal, 0, -1, Cqp, 0, 0)
01806                 );
01807                 
01808                 grpAppldChrgVec.push_back(
01809                     ComputeQMPntChrg(chrgPos2, trgChrg2, trgIndxLocal, 0, -1, Cqm, 0, 0)
01810                 );
01811                 
01812                 } break;
01813                 
01814                 // Redistributed Charge and Dipole scheme.
01815                 case QMSCHEMERCD:
01816                 {
01817                 
01818 //                 grpAppldChrgVec[trgIndxLocal].charge -= mmCharge ;
01819                 
01820                 // Redistributed Charge and Dipole method (DOI: 10.1021/ct100530r)
01821                 // Here we create a dipole to counter act the charge movement
01822                 // we created by moving parts of the MM charge to target MM atoms.
01823                 // 
01824                 // The MM1 charge is equally divided and subtracted from all MM2 atoms.
01825                 // One point charge is created in the midpoint of the MM1-MM2 bond.
01826                 // 
01827                 // Below we see the diagram of the positions of new charges along
01828                 // the direction of the bond between the MM1 atom and a target 
01829                 // MM2 atom, wiht respect to the new Dummy atom (a Hydrogen).
01830                 // 
01831                 //  QM --------  MM1(p0) -------------- MM2
01832                 //  QM ------ H ------------ (+2*p0) - (MM2 -p0)
01833                 
01834                 Position trgPos = grpAppldChrgVec[trgIndxLocal].position ;
01835                 
01836                 // We now store in MM2 to which PC it is bound (the index for MM1).
01837                 // This will be used to decompose virtual PC forces onto MM1 and MM2.
01838                 grpAppldChrgVec[trgIndxLocal].qmGrpID = -2;
01839                 grpAppldChrgVec[trgIndxLocal].homeIndx = mmIndex;
01840                 
01841                 // We create a new point charge at the same position so that
01842                 // the fraction of charge from the MM1 atom is "added" to the
01843                 // MM2 position, but without actually changing the charge of 
01844                 // the original MM2 atom. This trick helps keeping the original 
01845                 // charge for electrostatic calculations after the new charges 
01846                 // of QM atoms is received from the QM software.
01847                 Real Cq0 = 1.0;
01848                 grpAppldChrgVec.push_back(
01849                     // Cq is stored in the distance slot, since it has not been computed
01850                     // for a virtual charge.
01851                     ComputeQMPntChrg(trgPos, -1*mmCharge, trgIndxLocal, 0, -1, Cq0, 0, 0)
01852                 );
01853                 
01854                 Vector bondVec = trgPos - mmPos ;
01855                 
01856                 // For Cq-plus virtual charge
01857                 Real Cq1 = 0.5;
01858                 Vector bondVec1 = bondVec*Cq1 ;
01859                 
01860                 Position chrgPos1 = mmPos + bondVec1;
01861                 
01862                 BigReal trgChrg1 = 2*mmCharge;
01863                 
01864                 grpAppldChrgVec.push_back(
01865                     ComputeQMPntChrg(chrgPos1, trgChrg1, trgIndxLocal, 0, -1, Cq1, 0, 0)
01866                 );
01867                 
01868                 
01869                 
01870                 } break;
01871                 
01872                 // The following schemes simply make the charges of the 
01873                 // MM atom(s) equal to zero. No redistribution or 
01874                 // artificial dipole is attmepted.
01875                 //
01876                 // For Z1, only the MM1 atom has its charge set to zero.
01877                 // For Z2, the MM1 and all MM2 atom charges are changed.
01878                 // For Z3, the MM1, all MM2 and all MM3 atom charges are changed.
01879                 //
01880                 // All modification are local ONLY. They do not change the atom's
01881                 // charge for the classical side of the simulation. Only the QM
01882                 // side ceases to see the electrostatic influence of these atoms,
01883                 // and they cease to see the QM region electrostatics as well.
01884                 
01885                 // Z1, Z2 and Z3 schemes do the same thing. The only difference 
01886                 // is the target list, which is created in MoleculeQM.C.
01887                 case QMSCHEMEZ1:
01888                 
01889                 case QMSCHEMEZ2:
01890                 
01891                 case QMSCHEMEZ3:
01892                     grpAppldChrgVec[trgIndxLocal].charge = 0.0;
01893                     break;
01894             }
01895             
01896         }
01897         
01898         // We keep this "point charge" so we can calculate forces on it later
01899         // but it will not influence the QM system.
01900         // We use the qmGrpID variable to send the message that this point charge 
01901         // should be ignored since this variable will not be relevant anymore.
01902         // All point charges gathered here are for a specific qmGroup.
01903         grpAppldChrgVec[mmIndex].qmGrpID = -1 ;
01904     }
01905     
01906     return ;
01907 }
01908 
01909 
01910 void ComputeQMMgr::recvPntChrg(QMPntChrgMsg *msg) {
01911     
01912     // All the preparation that used to be in this function was moved to 
01913     // recvPartQM, which is called first in rank zero.
01914     
01915     if (noPC) {
01916         // Even if we have QM-MM bonds, the point charge messages 
01917         // are handled in recvPartQM
01918         delete msg;
01919     }
01920     else {
01921         pntChrgCoordMsgs[numArrivedPntChrgMsg] = msg;
01922         ++numArrivedPntChrgMsg;
01923 
01924         if ( numArrivedPntChrgMsg < numSources ) 
01925             return;
01926     }
01927     
01928     // Resets counters for next step.
01929     numRecQMRes = 0;
01930     
01931     totalEnergy = 0;
01932     
01933     for ( int k=0; k<3; ++k )
01934         for ( int l=0; l<3; ++l )
01935             totVirial[k][l] = 0;
01936     
01937     // ALL DATA ARRIVED --- CALCULATE FORCES
01938     
01939     const char *const *const elementArray = molPtr->get_qmElements() ;
01940     const char *const *const dummyElmArray = molPtr->get_qmDummyElement();
01941     const int *const qmGrpSize = molPtr->get_qmGrpSizes();
01942     
01943     const BigReal *const dummyBondVal = molPtr->get_qmDummyBondVal();
01944     const int *const grpNumBonds = molPtr->get_qmGrpNumBonds() ;
01945     const int *const *const qmMMBond = molPtr->get_qmMMBond() ;
01946     const int *const *const qmGrpBonds = molPtr->get_qmGrpBonds() ;
01947     const int *const *const qmMMBondedIndx = molPtr->get_qmMMBondedIndx() ;
01948     const int *const *const chargeTarget = molPtr->get_qmMMChargeTarget() ;
01949     const int *const numTargs = molPtr->get_qmMMNumTargs() ;
01950     
01951     BigReal constants = COULOMB * simParams->nonbondedScaling / (simParams->dielectric) ;
01952     // COULOMB is in kcal*Angs/(mol*e^2)
01953 //     BigReal constants = COULOMB ;
01954     
01955     if ( qmPCFreq > 0 ) {
01956         DebugM(4,"Using point charge stride of " << qmPCFreq << "\n")
01957         // In case we are skiping steps between PC selection, store a main list
01958         // in rank zero for future steps. Then we only update positions untill
01959         // we reach a new "qmPCFreq" step, when a new PC selection is made.
01960         
01961         if (timeStep % qmPCFreq == 0) {
01962             DebugM(4,"Loading a new selection of PCs.\n")
01963             
01964             // We only re-set this arrya in a step where a new PC selection is made.
01965             pntChrgMsgVec.clear();
01966             for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
01967                 // Accumulates the message point charges into a local vector.
01968                 for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
01969                     pntChrgMsgVec.push_back(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
01970                 }
01971             }
01972             
01973             // This fast array is created to send the entire PC IDs list to all
01974             // PEs in the next step.
01975             pcIDListSize = pntChrgMsgVec.size();
01976             pcIDList = new int[pcIDListSize] ;
01977             for (size_t i=0; i<pcIDListSize; i++) {
01978                 pcIDList[i] = pntChrgMsgVec[i].id;
01979                 
01980                 // Loads home indexes of MM atoms received as point charges.
01981                 // Each force receives the home index of its atom with respect to the 
01982                 // local set of atoms in each node.
01983                 force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
01984             }
01985         }
01986         else {
01987             DebugM(4,"Updating position/homeIdex of old PC selection.\n")
01988             
01989             // We build a sorted array so that we can quickly access the unordered
01990             // data we just received, and update positions of the same selection
01991             // of point charges.
01992             pcDataSort.clear();
01993             for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
01994                 // Accumulates the message point charges into a local sorted array.
01995                 for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
01996                     pcDataSort.load(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
01997                 }
01998             }
01999             pcDataSort.sort();
02000             
02001             iout << "Loaded new positions in sorted array: " << pcDataSort.size() << "\n" << endi;
02002             
02003             // If we are using a set of point charges that was selected in a
02004             // previous step, we update the PC POSITION ONLY.
02005             for (size_t i=0; i<pntChrgMsgVec.size(); i++) {
02006                 
02007                 auto pntr = pcDataSort.find(pntChrgMsgVec[i]);
02008                 
02009                 pntChrgMsgVec[i].position = pntr->position ;
02010                 pntChrgMsgVec[i].homeIndx = pntr->homeIndx ;
02011                 
02012                 // Loads home indexes of MM atoms received as point charges.
02013                 // Each force receives the home index of its atom with respect to the 
02014                 // local set of atoms in each node.
02015                 force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
02016             }
02017         }
02018     }
02019     else {
02020         DebugM(4,"Updating PCs at every step.\n")
02021         
02022         pntChrgMsgVec.clear();
02023         for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
02024             // Accumulates the message point charges into a local vector.
02025             for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
02026                 pntChrgMsgVec.push_back(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
02027             }
02028         }
02029         
02030         // Loads home indexes of MM atoms received as point charges.
02031         for (size_t i=0; i<pntChrgMsgVec.size(); i++) {
02032             // Each force receives the home index of its atom with respect to the 
02033             // local set of atoms in each node.
02034             
02035             #ifdef DEBUG_QM
02036             if (force[pntChrgMsgVec[i].id].homeIndx != -1 
02037                 and force[pntChrgMsgVec[i].id].homeIndx != pntChrgMsgVec[i].homeIndx
02038             ) {
02039                 iout << iERROR << "Overloading point charge " 
02040                 << pntChrgMsgVec[i].id << "; home index: " 
02041                 << force[pntChrgMsgVec[i].id].homeIndx << " with " << pntChrgMsgVec[i].homeIndx
02042                 << "\n" << endi ;
02043                 NAMD_die("Error preparing QM simulations.");
02044             }
02045             #endif
02046             
02047             force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
02048         }
02049     }
02050     
02051     // Reset counter for next timestep
02052     numArrivedPntChrgMsg = 0;
02053     
02054     DebugM(4,"A total of " << pntChrgMsgVec.size() 
02055     << " point charges arrived." << std::endl);
02056     
02057     DebugM(4,"Starting QM groups processing.\n");
02058     
02059     QMAtmVec grpQMAtmVec;
02060     QMPCVec grpPntChrgVec;
02061     
02062     // Final set of charges, created or not, that is sent to the QM software.
02063     // This set will depend on how QM-MM bonds are processed and presented to the
02064     // QM region.
02065     QMPCVec grpAppldChrgVec;
02066     
02067     // Vector of dummy atoms created to treat QM-MM bonds.
02068     std::vector<dummyData> dummyAtoms ;
02069 
02070     // Initializes the loop for receiving the QM results.
02071     thisProxy[0].recvQMResLoop() ;
02072     
02073     // Iterator for target PE where QM simulations will run.
02074     int peIter = 0;
02075     
02076     for (int grpIter = 0; grpIter < numQMGrps; grpIter++) {
02077         
02078         grpQMAtmVec.clear();
02079         grpPntChrgVec.clear();
02080         grpAppldChrgVec.clear();
02081         dummyAtoms.clear();
02082         
02083         DebugM(4,"Calculating QM group " << grpIter +1 
02084         << " (ID: " << grpID[grpIter] << ")." << std::endl);
02085         
02086         DebugM(4,"Compiling Config Lines into one string for message...\n");
02087         
02088         // This will hold a big sting with all configuration lines the user supplied.
02089         int lineIter = 0 ;
02090         std::string configLines ;
02091         StringList *current = Node::Object()->configList->find("QMConfigLine");
02092         for ( ; current; current = current->next ) {
02093             std::string confLineStr(current->data);
02094             
02095             // In case we need to add charges to MOPAC command line.
02096             if (simParams->qmFormat == QMFormatMOPAC && simParams->qmMOPACAddConfigChrg && lineIter == 0) {
02097                 std::ostringstream itosConv ;
02098                 itosConv << grpChrg[grpIter] ;
02099                 confLineStr.append( " CHARGE=" );
02100                 confLineStr.append( itosConv.str() );
02101                 
02102             }
02103             
02104             configLines.append(confLineStr);
02105             configLines.append("\n");
02106             
02107             lineIter++;
02108         }
02109         
02110         DebugM(4,"Determining point charges...\n");
02111         
02112         Real qmTotalCharge = 0;
02113         // Loads the QM atoms for this group.
02114         for (int qmIt=0; qmIt<numQMAtms; qmIt++){
02115             if (qmCoord[qmIt].qmGrpID == grpID[grpIter]) {
02116                 grpQMAtmVec.push_back(qmCoord[qmIt]);
02117                 qmTotalCharge += qmCoord[qmIt].charge;
02118             }
02119         }
02120         if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
02121             qmTotalCharge = roundf(qmTotalCharge) ;
02122         }
02123         
02124         // Sorts the vector so that the QM message is always built with the same order of atoms.
02125         std::sort(grpQMAtmVec.begin(), grpQMAtmVec.end(), custom_ComputeQMAtom_Less);
02126         
02127         Real pcTotalCharge = 0;
02128         // Loads the point charges to a local vector for this QM group.
02129         for (auto pntChrgIt = pntChrgMsgVec.begin(); 
02130              pntChrgIt != pntChrgMsgVec.end(); pntChrgIt++) {
02131             if ((*pntChrgIt).qmGrpID == grpID[grpIter] ) {
02132                 grpPntChrgVec.push_back( (*pntChrgIt) );
02133                 pcTotalCharge += (*pntChrgIt).charge;
02134             }
02135         }
02136         if ((fabsf(roundf(pcTotalCharge) - pcTotalCharge) <= 0.001f) ) {
02137             pcTotalCharge = roundf(pcTotalCharge) ;
02138         }
02139         
02140         #ifdef DEBUG_QM
02141         if (grpQMAtmVec.size() != qmGrpSize[grpIter]) {
02142             iout << iERROR << "The number of QM atoms received for group " << grpID[grpIter]
02143             << " does not match the expected: " << grpQMAtmVec.size()
02144             << " vs " << qmGrpSize[grpIter] << "\n" << endi ;
02145             
02146             NAMD_die("Error processing QM group.");
02147         }
02148         #endif
02149         
02150         DebugM(1,"Found " << grpPntChrgVec.size() << " point charges for QM group " 
02151         << grpIter << " (ID: " << grpID[grpIter] << "; Num QM atoms: " 
02152         << grpQMAtmVec.size() <<  "; Num QM-MM bonds: " 
02153         << grpNumBonds[grpIter] << ")" << std::endl);
02154         
02155         DebugM(1,"Total QM charge: " << qmTotalCharge 
02156         << "; Total point-charge charge: " << pcTotalCharge << std::endl);
02157         
02158         // If we have a frequency for LSS update, check if we shoudl do it in 
02159         // the current time step.
02160         if ( qmLSSFreq > 0 && ((timeStep + 1) % qmLSSFreq == 0 )) {
02161             lssUpdate(grpIter, grpQMAtmVec, grpPntChrgVec);
02162         }
02163         
02164         // This function checks data and treats the charge (and existence) of
02165         // the MM atoms in and around QM-MM bonds. It is only executed in 
02166         // electrostatic embeding QM/MM simulations.
02167         if (! noPC ) {
02168             procBonds(grpNumBonds[grpIter], qmGrpBonds[grpIter], 
02169                      qmMMBondedIndx[grpIter], 
02170                      chargeTarget, numTargs, 
02171                      grpPntChrgVec, grpAppldChrgVec) ;
02172         }
02173         else {
02174             grpAppldChrgVec = grpPntChrgVec;
02175         }
02176         
02177         // For future use, we get the pairs of indexes of QM atoms and MM atoms which are
02178         // bound in QM-MM bonds.
02179         std::vector< std::pair<int,int> > qmPCLocalIndxPairs ;
02180         
02181         // Create and position dummy atoms.
02182         Position mmPos(0), qmPos(0);
02183         for (int dummyIt = 0; dummyIt < grpNumBonds[grpIter]; dummyIt++) {
02184             
02185             int qmMMBondIndx = qmGrpBonds[grpIter][dummyIt] ;
02186             
02187             BigReal bondVal = dummyBondVal[qmMMBondIndx] ;
02188             
02189             int mmAtmIndx = qmMMBond[qmMMBondIndx][0] ;
02190             int qmAtmIndx = qmMMBond[qmMMBondIndx][1] ;
02191             
02192             // Sicne we don't know in which patch/node the QM atom is, or the 
02193             // order in which they will arrive in rank zero, we have
02194             // no direct index to it.
02195             #ifdef DEBUG_QM
02196             bool missingQM = true, missingMM = true;
02197             #endif
02198             size_t qmIt ;
02199             for (qmIt=0; qmIt<grpQMAtmVec.size(); qmIt++){
02200                 if (grpQMAtmVec[qmIt].id == qmAtmIndx) {
02201                     qmPos = grpQMAtmVec[qmIt].position;
02202                     
02203                     #ifdef DEBUG_QM
02204                     missingQM = false;
02205                     #endif
02206                     
02207                     break;
02208                 }
02209             }
02210             // The same is true about the MM atom to which the QM atom is bound,
02211             // we must look
02212             size_t pcIt;
02213             for (pcIt=0; pcIt < grpPntChrgVec.size(); pcIt++) {
02214                 if (grpPntChrgVec[pcIt].id == mmAtmIndx) {
02215                     mmPos = grpPntChrgVec[pcIt].position ;
02216                     
02217                     #ifdef DEBUG_QM
02218                     missingMM = false;
02219                     #endif
02220                     
02221                     break;
02222                 }
02223             }
02224             
02225             qmPCLocalIndxPairs.push_back(std::pair<int,int>(qmIt, pcIt) );
02226             
02227             #ifdef DEBUG_QM
02228             // Checks if the MM atom was included as a 
02229             // point charge due proximity to a QM atom, and if the QM atom arrived.
02230             if ( missingMM or missingQM ) {
02231                 // If it wasn't, there is an error somwhere.
02232                 
02233                 if (missingMM) {
02234                     iout << iERROR << "The MM atom " << mmAtmIndx
02235                     << " is bound to a QM atom, but it was not selected as a poitn charge."
02236                     << " Check your cutoff radius!\n" << endi ;
02237                     
02238                     NAMD_die("Error in QM-MM bond processing.");
02239                 }
02240                 if (missingQM) {
02241                     iout << iERROR << "The QM atom " << qmAtmIndx
02242                     << " is bound to an MM atom, but it was not sent to rank zero for processing."
02243                     << " Check your configuration!\n" << endi ;
02244                     
02245                     NAMD_die("Error in QM-MM bond processing.");
02246                 }
02247             }
02248             #endif
02249             
02250             Vector bondVec = mmPos - qmPos ;
02251             
02252             if (bondValType == QMLENTYPE) {
02253                 // If a length is defined by the user, or a default len
02254                 // is used, we calculate the unit vector for the displacement
02255                 // and multiply by the desired length in order 
02256                 // to get the final dummy atom position relative to the
02257                 // QM atom.
02258                 bondVec = bondVec.unit() ;
02259                 bondVec *= bondVal ;
02260             }
02261             else if (bondValType == QMRATIOTYPE) {
02262                 // If distance a ratio was defined by the user, then
02263                 // the displacement vector is multiplied by that ratio
02264                 // to get the final dummy atom position relative to the
02265                 // QM atom.
02266                 bondVec *= bondVal ;
02267             }
02268             
02269             Position dummyPos = qmPos + bondVec;
02270             
02271             DebugM(1,"Creating dummy atom " << dummyPos << " ; QM ind: " 
02272             << qmIt << " ; PC ind: " << pcIt << std::endl);
02273             
02274             dummyAtoms.push_back(dummyData(dummyPos, qmIt, qmMMBondIndx)) ;
02275             
02276         }
02277         
02278         DebugM(3, "Creating data for " << grpQMAtmVec.size() << " QM atoms " 
02279         << dummyAtoms.size() << " dummy atoms " << grpPntChrgVec.size()
02280         << " real point charges and " << grpAppldChrgVec.size() - grpPntChrgVec.size()
02281         << " virtual point charges\n") ;
02282         
02283         int dataSize = grpQMAtmVec.size() + dummyAtoms.size() + grpAppldChrgVec.size();
02284         QMGrpCalcMsg *msg = new (dataSize, configLines.size(), 0) QMGrpCalcMsg;
02285         msg->timestep = timeStep;
02286         msg->grpIndx = grpIter;
02287         msg->peIter = peIter;
02288         msg->charge = grpChrg[grpIter];
02289         msg->multiplicity = grpMult[grpIter];
02290         msg->numQMAtoms = grpQMAtmVec.size();
02291         msg->numAllAtoms = grpQMAtmVec.size() + dummyAtoms.size();
02292         msg->numRealPntChrgs = grpPntChrgVec.size(); // The original set of classical atoms.
02293         msg->numAllPntChrgs = grpAppldChrgVec.size(); // The extended set with virtual point charges.
02294         msg->secProcOn = simParams->qmSecProcOn ;
02295         msg->constants = constants;
02296         msg->PMEOn = simParams->PMEOn ;
02297         if (msg->PMEOn)
02298             msg->PMEEwaldCoefficient = simParams->PMEEwaldCoefficient ;
02299         msg->switching = simParams->qmPCSwitchOn;
02300         msg->switchType = simParams->qmPCSwitchType;
02301         msg->cutoff = simParams->cutoff;
02302         msg->swdist = simParams->switchingDist;
02303         msg->pcScheme = simParams->qmPCScheme;
02304         msg->qmAtmChrgMode = simParams->qmChrgMode;
02305         
02306         strncpy(msg->baseDir, simParams->qmBaseDir, 256);
02307         strncpy(msg->execPath, simParams->qmExecPath, 256);
02308         if (msg->secProcOn)
02309             strncpy(msg->secProc, simParams->qmSecProc, 256);
02310         
02311         if (simParams->qmPrepProcOn && (timeStep == simParams->firstTimestep)) {
02312             msg->prepProcOn = true;
02313             strncpy(msg->prepProc, simParams->qmPrepProc, 256);
02314         } else
02315             msg->prepProcOn = false;
02316         
02317         QMAtomData *dataP = msg->data;
02318         
02319         for (int i=0; i<grpQMAtmVec.size(); i++) {
02320             dataP->position = grpQMAtmVec[i].position ;
02321             dataP->charge = grpQMAtmVec[i].charge ;
02322             dataP->id = grpQMAtmVec[i].id ;
02323             dataP->bountToIndx = -1;
02324             dataP->type = QMATOMTYPE_QM ;
02325             strncpy(dataP->element,elementArray[grpQMAtmVec[i].id],3);
02326             dataP++;
02327         }
02328         
02329         for (int i=0; i<dummyAtoms.size(); i++) {
02330             dataP->position = dummyAtoms[i].pos ;
02331             dataP->charge = 0 ;
02332             dataP->id = -1 ;
02333             dataP->bountToIndx = dummyAtoms[i].qmInd ;
02334             dataP->type = QMATOMTYPE_DUMMY ;
02335             strncpy(dataP->element,dummyElmArray[dummyAtoms[i].bondIndx],3);
02336             dataP++;
02337         }
02338         
02339         for (int i=0; i<grpAppldChrgVec.size(); i++) {
02340             dataP->position = grpAppldChrgVec[i].position ;
02341             dataP->charge = grpAppldChrgVec[i].charge ;
02342             // Point charges created to handle QM-MM bonds will have an id of -1.
02343             dataP->id = grpAppldChrgVec[i].id ;
02344             dataP->bountToIndx = -1;
02345             dataP->dist = grpAppldChrgVec[i].dist ;
02346             // If we are loading the classical atoms' charges
02347             // the point charge type is 1, unless it is from an 
02348             // atom which is bound to a QM atom.
02349             if (i < grpPntChrgVec.size()) {
02350                 if (grpAppldChrgVec[i].qmGrpID == -1) {
02351                     dataP->type = QMPCTYPE_IGNORE ;
02352                 }
02353                 else if (grpAppldChrgVec[i].qmGrpID == -2) {
02354                     dataP->type = QMPCTYPE_CLASSICAL ;
02355                     dataP->bountToIndx = grpAppldChrgVec[i].homeIndx;
02356                 }
02357                 else {
02358                     dataP->type = QMPCTYPE_CLASSICAL ;
02359                 }
02360             }
02361             else {
02362                 // Extra charges are created to handle QM-MM bonds (if they exist).
02363                 dataP->type = QMPCTYPE_EXTRA ;
02364                 dataP->bountToIndx = grpAppldChrgVec[i].id;
02365             }
02366             dataP++;
02367         }
02368         
02369         QMAtomData *qmP = msg->data ;
02370         QMAtomData *pcP = msg->data + msg->numAllAtoms ;
02371         
02372         // With this, every QM atom knows to which MM atom is is bound, 
02373         // and vice versa. This will be usefull later on to prevent them from
02374         // feeling eachother's electrostatic charges AND to place the dummy
02375         // atom forces on the "real" atoms that form the bond.
02376         for( auto vecPtr  = qmPCLocalIndxPairs.begin(); 
02377                   vecPtr != qmPCLocalIndxPairs.end(); 
02378                   vecPtr++ ) {
02379             
02380             int qmLocInd = (*vecPtr).first;
02381             int pcLocInd = (*vecPtr).second;
02382             
02383             qmP[qmLocInd].bountToIndx = pcLocInd ;
02384             pcP[pcLocInd].bountToIndx = qmLocInd ;
02385         }
02386         
02387         
02388         strcpy(msg->configLines, configLines.c_str());
02389         
02390         if (cSMDon)
02391             calcCSMD(grpIter, msg->numQMAtoms, qmP, cSMDForces) ;
02392                 
02393         int targetPE = qmPEs[peIter] ;
02394         
02395         DebugM(4,"Sending QM group " << grpIter << " (ID " << grpID[grpIter] 
02396         << ") to PE " << targetPE << std::endl);
02397         
02398         switch (simParams->qmFormat) {
02399             // Creates the command line that will be executed according to the 
02400             // chosen QM software, as well as the input file with coordinates.
02401             case QMFormatORCA:
02402                 QMProxy[targetPE].calcORCA(msg) ;
02403                 break;
02404             
02405             case QMFormatMOPAC:
02406                 QMProxy[targetPE].calcMOPAC(msg) ;
02407                 break;
02408             
02409             case QMFormatUSR:
02410                 QMProxy[targetPE].calcUSR(msg) ;
02411                 break;
02412         }
02413         
02414         peIter++;
02415         
02416         if (peIter == qmPEs.size())
02417             peIter = 0;
02418     }
02419     
02420 }
02421 
02422 void ComputeQMMgr::storeQMRes(QMGrpResMsg *resMsg) {
02423     
02424 //     iout << iINFO << "Storing QM results for region " << resMsg->grpIndx  
02425 //     << " (ID: "  << grpID[resMsg->grpIndx] 
02426 //     << ") with original energy: " << endi;
02427 //     std::cout << std::fixed << std::setprecision(6) << resMsg->energyOrig << endi;
02428 //     iout << " Kcal/mol\n" << endi;
02429     
02430 //     if (resMsg->energyCorr != resMsg->energyOrig) {
02431 //         iout << iINFO << "PME corrected energy: " << endi;
02432 //         std::cout << std::fixed << std::setprecision(6) << resMsg->energyCorr << endi;
02433 //         iout << " Kcal/mol\n" << endi;
02434 //     }
02435     
02436     if ( (timeStep % simParams->qmEnergyOutFreq ) == 0 ) {
02437         iout << QMETITLE(timeStep);
02438         iout << FORMAT(grpID[resMsg->grpIndx]);
02439         iout << FORMAT(resMsg->energyOrig);
02440         if (resMsg->energyCorr != resMsg->energyOrig) iout << FORMAT(resMsg->energyCorr);
02441         iout << "\n\n" << endi;
02442     }
02443     
02444     totalEnergy += resMsg->energyCorr ;
02445     
02446     for ( int k=0; k<3; ++k )
02447         for ( int l=0; l<3; ++l )
02448             totVirial[k][l] += resMsg->virial[k][l];
02449     
02450     QMForce *fres = resMsg->force ;
02451     Real qmTotalCharge = 0;
02452     
02453     for (int i=0; i<resMsg->numForces; i++) {
02454         
02455         force[ fres[i].id ].force += fres[i].force;
02456         
02457         // Indicates the result is a QM atom, and we should get its updated charge.
02458         if (fres[i].replace == 1) {
02459             force[ fres[i].id ].charge =  fres[i].charge;
02460             qmTotalCharge += fres[i].charge;
02461         }
02462     }
02463     
02464     if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
02465         qmTotalCharge = roundf(qmTotalCharge) ;
02466     }
02467     
02468     // In case we are calculating cSMD forces, apply them now.
02469     if (cSMDon) {
02470         for ( int i=0; i < cSMDindxLen[resMsg->grpIndx]; i++ ) {
02471             int cSMDit = cSMDindex[resMsg->grpIndx][i];
02472             force[ cSMDpairs[cSMDit][0] ].force += cSMDForces[cSMDit];
02473         }
02474     }
02475 
02476     DebugM(4,"QM total charge received is " << qmTotalCharge << std::endl);
02477     
02478     DebugM(4,"Current accumulated energy is " << totalEnergy << std::endl);
02479     
02480     numRecQMRes++;
02481     
02482     delete resMsg;
02483 }
02484 
02485 void ComputeQMMgr::procQMRes() {
02486     
02487     // Writes a DCD file with the charges of all QM atoms at a frequency 
02488     // defined by the user in qmOutFreq.
02489     if ( simParams->qmOutFreq > 0 && 
02490          timeStep % simParams->qmOutFreq == 0 ) {
02491         
02492         iout << iINFO << "Writing QM charge output at step " 
02493         << timeStep <<  "\n" << endi;
02494     
02495         Real *x = outputData, 
02496              *y = outputData + molPtr->get_numQMAtoms(), 
02497              *z = outputData + 2*molPtr->get_numQMAtoms();
02498         
02499         for (int qmIt=0; qmIt<numQMAtms; qmIt++){
02500             x[qmIt] = qmCoord[qmIt].id;
02501             y[qmIt] = force[qmCoord[qmIt].id].charge ;
02502             z[qmIt] = 0;
02503         }
02504         
02505         write_dcdstep(dcdOutFile, numQMAtms, x, y, z, 0) ;
02506     }
02507     
02508     // Writes a DCD file with the positions of all QM atoms at a frequency 
02509     // defined by the user in qmPosOutFreq.
02510     if ( simParams->qmPosOutFreq > 0 && 
02511          timeStep % simParams->qmPosOutFreq == 0 ) {
02512         
02513         iout << iINFO << "Writing QM position output at step " 
02514         << timeStep <<  "\n" << endi;
02515         
02516         SortedArray<idIndxStr> idIndx;
02517         
02518         for(int i=0; i<numQMAtms;i++) {
02519             idIndx.insert( idIndxStr(qmCoord[i].id, i) );
02520         }
02521         idIndx.sort();
02522         
02523         Real *x = outputData, 
02524              *y = outputData + molPtr->get_numQMAtoms(), 
02525              *z = outputData + 2*molPtr->get_numQMAtoms();
02526         
02527         for (int qmIt=0; qmIt<numQMAtms; qmIt++){
02528             x[qmIt] = qmCoord[idIndx[qmIt].indx].position.x;
02529             y[qmIt] = qmCoord[idIndx[qmIt].indx].position.y;
02530             z[qmIt] = qmCoord[idIndx[qmIt].indx].position.z;
02531         }
02532         
02533         write_dcdstep(dcdPosOutFile, numQMAtms, x, y, z, 0) ;
02534     }
02535     
02536     // distribute forces
02537     DebugM(4,"Distributing QM forces for all ranks.\n");
02538     for ( int j=0; j < numSources; ++j ) {
02539         
02540         DebugM(1,"Sending forces and charges to source " << j << std::endl);
02541         
02542         QMCoordMsg *qmmsg = 0;
02543         QMPntChrgMsg *pcmsg = 0 ;
02544         
02545         int totForces = 0;
02546         int sourceNode = -1;
02547         
02548         if (pntChrgCoordMsgs == NULL) {
02549             
02550             qmmsg = qmCoordMsgs[j];
02551             qmCoordMsgs[j] = 0;
02552             
02553             totForces = qmmsg->numAtoms ;
02554             
02555             sourceNode = qmmsg->sourceNode;
02556         }
02557         else {
02558             pcmsg = pntChrgCoordMsgs[j];
02559             pntChrgCoordMsgs[j] = 0;
02560             
02561             sourceNode = pcmsg->sourceNode;
02562             
02563             // Since we receive two different messages from nodes, there is no 
02564             // guarantee the two sets of messages will come in the same order.
02565             // Therefore, we match the messages by comaring their sourceNodes.
02566             for (int aux=0; aux<numSources; aux++) {
02567                 
02568                 if (qmCoordMsgs[aux] == 0)
02569                     continue;
02570                 
02571                 qmmsg = qmCoordMsgs[aux];
02572                 
02573                 if (qmmsg->sourceNode == sourceNode) {
02574                     qmCoordMsgs[aux] = 0;
02575                     break;
02576                 }
02577             }
02578             
02579             DebugM(1,"Building force mesage for rank " 
02580             << pcmsg->sourceNode << std::endl);
02581             
02582             totForces = qmmsg->numAtoms + pcmsg->numAtoms;
02583         }
02584         
02585         QMForceMsg *fmsg = new (totForces, subsArray.size(), 0) QMForceMsg;
02586         fmsg->PMEOn = simParams->PMEOn;
02587         fmsg->numForces = totForces;
02588         fmsg->numLssDat = subsArray.size();
02589         
02590         DebugM(1,"Loading QM forces.\n");
02591         
02592         // This iterator is used in BOTH loops!
02593         int forceIter = 0;
02594         
02595         for ( int i=0; i < qmmsg->numAtoms; ++i ) {
02596             fmsg->force[forceIter] = force[qmmsg->coord[i].id];
02597             forceIter++;
02598         }
02599         
02600         delete qmmsg;
02601         
02602         if (pntChrgCoordMsgs != NULL) {
02603             DebugM(1,"Loading PntChrg forces.\n");
02604             
02605             for ( int i=0; i < pcmsg->numAtoms; ++i ) {
02606                 fmsg->force[forceIter] = force[pcmsg->coord[i].id];
02607                 forceIter++;
02608             }
02609             
02610             delete pcmsg;
02611         }
02612         
02613         DebugM(1,"A total of " << forceIter << " forces were loaded." << std::endl);
02614         
02615         for ( int i=0; i < subsArray.size(); ++i ) {
02616             fmsg->lssDat[i] = subsArray[i];
02617         }
02618         
02619         #ifdef DEBUG_QM
02620         if (subsArray.size() > 0)
02621             DebugM(3,"A total of " << subsArray.size() << " LSS data structures were loaded." << std::endl);
02622         #endif
02623         
02624         if ( ! j ) {
02625             fmsg->energy = totalEnergy;
02626             for ( int k=0; k<3; ++k )
02627                 for ( int l=0; l<3; ++l )
02628                     fmsg->virial[k][l] = totVirial[k][l];
02629         } else {
02630             fmsg->energy = 0;
02631             for ( int k=0; k<3; ++k )
02632                 for ( int l=0; l<3; ++l )
02633                     fmsg->virial[k][l] = 0;
02634         }
02635         
02636         DebugM(4,"Sending forces...\n");
02637         
02638         QMProxy[sourceNode].recvForce(fmsg);
02639         
02640     }
02641     
02642     DebugM(4,"All forces sent from node zero.\n");
02643 }
02644 
02645 void ComputeQMMgr::recvForce(QMForceMsg *fmsg) {
02646     
02647     if (CkMyPe()) {
02648         for (int i=0; i<fmsg->numLssDat; i++) {
02649             subsArray.add(fmsg->lssDat[i]) ;
02650         }
02651     }
02652     
02653     QMCompute->saveResults(fmsg);
02654 }
02655 
02656 void ComputeQM::saveResults(QMForceMsg *fmsg) {
02657     
02658     ResizeArrayIter<PatchElem> ap(patchList);
02659     
02660     bool callReplaceForces = false;
02661     
02662     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
02663     const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
02664     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
02665     Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
02666     
02667     int totalNumbAtoms = 0;
02668     for (ap = ap.begin(); ap != ap.end(); ap++) {
02669         totalNumbAtoms += (*ap).p->getNumAtoms();
02670     }
02671     
02672     // This is kept to be deleted in the next step so that the message can be 
02673     // used in "replaceForces" routine later on, in case there are forces to 
02674     // be replaced. The "replaceForces" function uses but DOES NOT free the pointer,
02675     // so we free the data from the previous iteration and allocate a new one for
02676     // the current iteration (remember that the number of atoms can change in a 
02677     // home patch between iterations).
02678     if (oldForces != 0)
02679         delete [] oldForces;
02680     oldForces = new ExtForce[totalNumbAtoms] ;
02681     
02682     for (int i=0; i < totalNumbAtoms; ++i) {
02683         oldForces[i].force = Force(0) ;
02684     }
02685     
02686     DebugM(1,"Force array has been created and zeroed in rank " 
02687     << CkMyPe() << std::endl);
02688     
02689     DebugM(1,"Preparing " << fmsg->numForces << " forces in rank " 
02690     << CkMyPe() << std::endl);
02691     
02692     QMForce *results_ptr = fmsg->force;
02693     // Iterates over the received forces and place them on the full array.
02694     for (int i=0; i < fmsg->numForces; ++i, ++results_ptr) {
02695         // For now we may have more than one item in the message acting on the same 
02696         // atom, such as an MM atom which is a point charge for two or more QM regions.
02697         
02698         oldForces[results_ptr->homeIndx].force += results_ptr->force;
02699         oldForces[results_ptr->homeIndx].replace = results_ptr->replace;
02700         
02701         if (results_ptr->replace == 1)
02702             callReplaceForces = true;
02703         
02704         // If the atom is in a QM group, update its charge to the local (this homePatch)
02705         // copy of the qmAtmChrg array.
02706         if (qmAtomGroup[results_ptr->id] > 0 && fmsg->PMEOn) {
02707             
02708             // Loops over all QM atoms (in all QM groups) comparing their global indices
02709             for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
02710                 
02711                 if (qmAtmIndx[qmIter] == results_ptr->id) {
02712                     qmAtmChrg[qmIter] = results_ptr->charge;
02713                     break;
02714                 }
02715                 
02716             }
02717             
02718         }
02719         
02720     }
02721     
02722     DebugM(1,"Placing forces on force boxes in rank " 
02723     << CkMyPe() << std::endl);
02724     
02725     // Places the received forces on the force array for each patch.
02726     int homeIndxIter = 0;
02727     for (ap = ap.begin(); ap != ap.end(); ap++) {
02728         Results *r = (*ap).forceBox->open();
02729         Force *f = r->f[Results::normal];
02730         const CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
02731         int localNumAtoms = (*ap).p->getNumAtoms();
02732         
02733         for(int i=0; i<localNumAtoms; ++i) {
02734             
02735             f[i] += oldForces[homeIndxIter].force;
02736             
02737             ++homeIndxIter;
02738         }
02739         
02740         if ( callReplaceForces )
02741             (*ap).p->replaceForces(oldForces);
02742         
02743         (*ap).forceBox->close(&r);
02744         
02745     }
02746     
02747     DebugM(1,"Forces placed on force boxes in rank " 
02748     << CkMyPe() << std::endl);
02749     
02750     if (fmsg->PMEOn) {
02751         
02752         DebugM(1,"PME ON! Accessing PMEmgr in rank " << CkMyPe() << std::endl);
02753         
02754         ComputePmeMgr *mgrP = CProxy_ComputePmeMgr::ckLocalBranch(
02755             CkpvAccess(BOCclass_group).computePmeMgr) ;
02756         
02757         DebugM(4, "Initiating ComputePme::doQMWork on rank " << CkMyPe() << " over " 
02758             << getComputes(mgrP).size() << " pmeComputes." << std::endl) ;
02759         
02760         for ( int i=0; i< getComputes(mgrP).size(); ++i ) {
02761     //         WorkDistrib::messageEnqueueWork(pmeComputes[i]);
02762             getComputes(mgrP)[i]->doQMWork();
02763         }
02764     }
02765     
02766     DebugM(1,"Submitting reduction in rank " << CkMyPe() << std::endl);
02767     
02768     reduction->item(REDUCTION_MISC_ENERGY) += fmsg->energy;
02769     reduction->item(REDUCTION_VIRIAL_NORMAL_XX) += fmsg->virial[0][0];
02770     reduction->item(REDUCTION_VIRIAL_NORMAL_XY) += fmsg->virial[0][1];
02771     reduction->item(REDUCTION_VIRIAL_NORMAL_XZ) += fmsg->virial[0][2];
02772     reduction->item(REDUCTION_VIRIAL_NORMAL_YX) += fmsg->virial[1][0];
02773     reduction->item(REDUCTION_VIRIAL_NORMAL_YY) += fmsg->virial[1][1];
02774     reduction->item(REDUCTION_VIRIAL_NORMAL_YZ) += fmsg->virial[1][2];
02775     reduction->item(REDUCTION_VIRIAL_NORMAL_ZX) += fmsg->virial[2][0];
02776     reduction->item(REDUCTION_VIRIAL_NORMAL_ZY) += fmsg->virial[2][1];
02777     reduction->item(REDUCTION_VIRIAL_NORMAL_ZZ) += fmsg->virial[2][2];
02778     reduction->submit();
02779     
02780     delete fmsg ;
02781 }
02782 
02783 void ComputeQMMgr::calcMOPAC(QMGrpCalcMsg *msg)
02784 {
02785     
02786     FILE *inputFile,*outputFile,*chrgFile;
02787     int iret;
02788     
02789     const size_t lineLen = 256;
02790     char *line = new char[lineLen];
02791     
02792     std::string qmCommand, inputFileName, outputFileName, pntChrgFileName;
02793     
02794     // For coulomb calculation
02795     BigReal constants = msg->constants;
02796     
02797     double gradient[3];
02798     
02799     DebugM(4,"Running MOPAC on PE " << CkMyPe() << std::endl);
02800     
02801     QMAtomData *atmP = msg->data ;
02802     QMAtomData *pcP = msg->data + msg->numAllAtoms ;
02803     
02804     // We create a pointer for PME correction, which may point to
02805     // a copy of the original point charge array, with unchanged charges,
02806     // or points to the original array in case no switching or charge
02807     // scheme is being used.
02808     QMAtomData *pcPpme = NULL;
02809     if (msg->switching) {
02810         
02811         if (msg->PMEOn)
02812             pcPpme = new QMAtomData[msg->numRealPntChrgs];
02813         
02814         pntChrgSwitching(msg, pcPpme) ;
02815     } else {
02816         pcPpme = pcP;
02817     }
02818     
02819     int retVal = 0;
02820     struct stat info;
02821     
02822     // For each QM group, create a subdirectory where files will be palced.
02823     std::string baseDir(msg->baseDir);
02824     std::ostringstream itosConv ;
02825     if ( CmiNumPartitions() > 1 ) {
02826         baseDir.append("/") ;
02827         itosConv << CmiMyPartition() ;
02828         baseDir += itosConv.str() ;
02829         itosConv.str("");
02830         itosConv.clear() ;
02831         
02832         if (stat(msg->baseDir, &info) != 0 ) {
02833             CkPrintf( "Node %d cannot access directory %s\n",
02834                       CkMyPe(), baseDir.c_str() );
02835             NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
02836         }
02837         else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
02838             DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
02839             retVal = mkdir(baseDir.c_str(), S_IRWXU);
02840         }
02841         
02842     }
02843     baseDir.append("/") ;
02844     itosConv << msg->peIter ;
02845     baseDir += itosConv.str() ;
02846     
02847     if (stat(msg->baseDir, &info) != 0 ) {
02848         CkPrintf( "Node %d cannot access directory %s\n",
02849                   CkMyPe(), baseDir.c_str() );
02850         NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
02851     }
02852     else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
02853         DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
02854         retVal = mkdir(baseDir.c_str(), S_IRWXU);
02855     }
02856     
02857     #ifdef DEBUG_QM
02858     Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
02859     #endif
02860     
02861     inputFileName.clear();
02862     inputFileName.append(baseDir.c_str()) ;
02863     inputFileName.append("/qmmm_") ;
02864     inputFileName += itosConv.str() ;
02865     inputFileName.append(".input") ;
02866     
02867     // Opens file for coordinate and parameter input
02868     inputFile = fopen(inputFileName.c_str(),"w");
02869     if ( ! inputFile ) {
02870         iout << iERROR << "Could not open input file for writing: " 
02871         << inputFileName << "\n" << endi ;
02872         NAMD_err("Error writing QM input file.");
02873     }
02874     
02875     // Builds the command that will be executed
02876     qmCommand.clear();
02877     qmCommand.append("cd ");
02878     qmCommand.append(baseDir);
02879     qmCommand.append(" ; ");
02880     qmCommand.append(msg->execPath) ;
02881     qmCommand.append(" ") ;
02882     qmCommand.append(inputFileName) ;
02883     qmCommand.append(" > /dev/null 2>&1") ;
02884     
02885     // Builds the file name where MOPAC will place the gradient
02886     // This is also relative to the input file name.
02887     outputFileName = inputFileName ;
02888     outputFileName.append(".aux") ;
02889     
02890     if (msg->numAllPntChrgs) {
02891         // Builds the file name where we will write the point charges.
02892         pntChrgFileName.clear();
02893         pntChrgFileName.append(baseDir.c_str()) ;
02894         pntChrgFileName.append("/mol.in") ;
02895         
02896         chrgFile = fopen(pntChrgFileName.c_str(),"w");
02897         if ( ! chrgFile ) {
02898             iout << iERROR << "Could not open charge file for writing: " 
02899             << pntChrgFileName << "\n" << endi ;
02900             NAMD_err("Error writing point charge file.");
02901         }
02902         
02903         iret = fprintf(chrgFile,"\n%d %d\n", 
02904                        msg->numQMAtoms, msg->numAllAtoms - msg->numQMAtoms);
02905         if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
02906     }
02907     
02908     // writes configuration lines to the input file.
02909     std::stringstream ss(msg->configLines) ;
02910     std::string confLineStr;
02911     while (std::getline(ss, confLineStr) ) {
02912         confLineStr.append("\n");
02913         iret = fprintf(inputFile,confLineStr.c_str());
02914         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
02915     }
02916     
02917     
02918     iret = fprintf(inputFile,"\n");
02919     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
02920     
02921     DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " 
02922     << inputFileName.c_str() << " and " << msg->numAllPntChrgs 
02923     << " point charges in file " << pntChrgFileName.c_str() << "\n");
02924     
02925     // write QM and dummy atom coordinates to input file and
02926     // MM electric field from MM point charges.
02927     for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
02928         
02929         double x = atmP->position.x;
02930         double y = atmP->position.y;
02931         double z = atmP->position.z;
02932         
02933         iret = fprintf(inputFile,"%s %f 1 %f 1 %f 1\n",
02934                        atmP->element,x,y,z);
02935         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
02936         
02937         if (msg->numAllPntChrgs) {
02938             BigReal phi = 0;
02939             
02940             // The Electrostatic Potential is calculated according to
02941             // the QMMM Keyword documentation for MOPAC
02942             // http://openmopac.net/manual/QMMM.html
02943             pcP = msg->data + msg->numAllAtoms ;
02944             for ( size_t j=0; j < msg->numAllPntChrgs; ++j, ++pcP ) {
02945                 
02946                 // In case of QM-MM bonds, the charge of the MM1 atom is ignored
02947                 if (pcP->type == QMPCTYPE_IGNORE)
02948                     continue;
02949                 
02950                 double charge = pcP->charge;
02951                 
02952                 double xMM = pcP->position.x;
02953                 double yMM = pcP->position.y;
02954                 double zMM = pcP->position.z;
02955                 
02956                 BigReal rij = sqrt((x-xMM)*(x-xMM) +
02957                      (y-yMM)*(y-yMM) +
02958                      (z-zMM)*(z-zMM) ) ;
02959                 
02960                 phi += charge/rij ;
02961             }
02962             
02963             // We use the same Coulomb constant used in the rest of NAMD
02964             // instead of the suggested rounded 332 suggested by MOPAC.
02965             phi = phi*constants ;
02966             
02967             iret = fprintf(chrgFile,"%s %f %f %f %f\n",
02968                            atmP->element,x,y,z, phi);
02969             if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
02970         }
02971     }
02972     
02973     DebugM(4,"Closing input and charge file\n");
02974     
02975     if (msg->numAllPntChrgs)
02976         fclose(chrgFile);
02977     
02978     fclose(inputFile);
02979     
02980     if (msg->prepProcOn) {
02981         
02982         std::string prepProc(msg->prepProc) ;
02983         prepProc.append(" ") ;
02984         prepProc.append(inputFileName) ;
02985         iret = system(prepProc.c_str());
02986         if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
02987         if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
02988     }
02989     
02990     // runs QM command
02991     DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
02992     iret = system(qmCommand.c_str());
02993     
02994     if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
02995     if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
02996 
02997     if (msg->secProcOn) {
02998         
02999         std::string secProc(msg->secProc) ;
03000         secProc.append(" ") ;
03001         secProc.append(inputFileName) ;
03002         itosConv.str("");
03003         itosConv.clear() ;
03004         itosConv << msg->timestep ;
03005         secProc.append(" ") ;
03006         secProc += itosConv.str() ;
03007         
03008         iret = system(secProc.c_str());
03009         if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
03010         if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
03011     }
03012     
03013     // remove coordinate file
03014 //     iret = remove(inputFileName);
03015 //     if ( iret ) { NAMD_err(0); }
03016 
03017     // remove coordinate file
03018 //     iret = remove(pntChrgFileName);
03019 //     if ( iret ) { NAMD_err(0); }
03020     
03021     // opens output file
03022     DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
03023     outputFile = fopen(outputFileName.c_str(),"r");
03024     if ( ! outputFile ) {
03025         iout << iERROR << "Could not find QM output file!\n" << endi;
03026         NAMD_err(0); 
03027     }
03028     
03029     // Resets the pointers.
03030     atmP = msg->data ;
03031     pcP = msg->data + msg->numAllAtoms ;
03032     
03033     // Allocates the return message.
03034     QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
03035     resMsg->grpIndx = msg->grpIndx;
03036     resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
03037     resMsg->energyOrig = 0;
03038     resMsg->energyCorr = 0;
03039     for ( int k=0; k<3; ++k )
03040         for ( int l=0; l<3; ++l )
03041             resMsg->virial[k][l] = 0;
03042     QMForce *resForce = resMsg->force ;
03043     
03044     // We split the data into two pointers so we can "skip" the dummy atoms
03045     // which have no representation in NAMD.
03046     for (int i=0; i<resMsg->numForces; i++) {
03047         resForce[i].force = 0;
03048         resForce[i].charge = 0 ;
03049         if (i < msg->numQMAtoms) {
03050             // We use the replace field to indicate QM atoms,
03051             // which will have charge information.
03052             resForce[i].replace = 1 ;
03053             resForce[i].id = atmP->id;
03054             atmP++;
03055         }
03056         else {
03057             // We use the replace field to indicate QM atoms,
03058             // which will have charge information.
03059             resForce[i].replace = 0 ;
03060             resForce[i].id = pcP->id;
03061             pcP++;
03062         }
03063     }
03064     
03065     // Resets the pointers.
03066     atmP = msg->data ;
03067     pcP = msg->data + msg->numAllAtoms ;
03068     
03069     // Reads the data form the output file created by the QM software.
03070     // Gradients over the QM atoms, and Charges for QM atoms will be read.
03071     
03072     size_t atmIndx = 0, gradCount = 0;
03073     Bool gradFields = false, chargeFields = false;
03074     Bool chargesRead = false, gradsRead = false;
03075     while ( fgets(line, lineLen, outputFile) != NULL){
03076         // We loop over the lines of the output file untill we find
03077         // the line that initiates the "atom charges" lines. Then
03078         // we read all lines and expect to get one or more charges 
03079         // per line, separated by space(s), untill we find a line that
03080         // initiates the "gradients", then once more, we expect several
03081         // numbers separated by space(s). When the "overlap matrix" 
03082         // string is found, we break the loop and stop reading the file.
03083         
03084         // if ( strstr(line,"TOTAL_ENERGY") != NULL ) {
03085         if ( strstr(line,"HEAT_OF_FORMATION") != NULL ) {
03086             
03087             char strEnergy[14], *endPtr ;
03088             
03089             //strncpy(strEnergy, line + 17, 13) ;
03090             strncpy(strEnergy, line + 28, 13) ;
03091             strEnergy[13] = '\0';
03092             
03093             // We have to convert the notation from FORTRAN double precision to
03094             // the natural, normal, reasonable and not terribly out dated format.
03095             resMsg->energyOrig = strtod(strEnergy, &endPtr);
03096             if ( *endPtr == 'D' ) {
03097                 *endPtr = 'e';
03098                 resMsg->energyOrig = strtod(strEnergy, &endPtr);
03099             }
03100             
03101             // In MOPAC, the total energy is given in EV, so we convert to Kcal/mol
03102 //             resMsg->energyOrig *= 23.060 ; // We now read Heat of Formation, which is given in Kcal/Mol
03103             
03104 //             DebugM(4,"Reading QM energy from file: " << resMsg->energyOrig << "\n");
03105             
03106             resMsg->energyCorr = resMsg->energyOrig;
03107             
03108             continue;
03109         }
03110         
03111         if ( strstr(line,"ATOM_CHARGES") != NULL ) {
03112             gradFields = false;
03113             chargeFields = true;
03114             atmIndx = 0;
03115             
03116             // If the user asked for charges NOT to be read from MOPAC,
03117             // we skip the charge reading step.
03118             if (msg->qmAtmChrgMode == QMCHRGNONE) {
03119                 chargeFields = false;
03120                 atmIndx = msg->numAllAtoms;
03121             }
03122             else {
03123                 chargeFields = true;
03124                 atmIndx = 0;
03125             }
03126             
03127             // Now we expect the following line(s) to have atom charges
03128             continue;
03129         }
03130         
03131         if ( strstr(line,"GRADIENTS") != NULL ) {
03132             
03133             // Now that all charges have been read, checks if the 
03134             // numbers match
03135             if (atmIndx != msg->numAllAtoms) {
03136                 NAMD_die("Error reading QM forces file. Wrong number of atom charges");
03137             }
03138             
03139             chargesRead = true;
03140             
03141             // Now we expect the following line(s) to have gradients
03142             chargeFields = false ;
03143             gradFields = true;
03144             gradCount = 0;
03145             atmIndx = 0;
03146             
03147             continue;
03148         }
03149         
03150         if ( strstr(line,"OVERLAP_MATRIX") != NULL ) {
03151             
03152             // Now that all gradients have been read, checks if the 
03153             // numbers match
03154             if (atmIndx != msg->numAllAtoms) {
03155                 NAMD_die("Error reading QM forces file. Wrong number of gradients");
03156             }
03157             
03158             gradsRead = true;
03159             
03160             // Nothing more to read from the ".aux" file
03161             break;
03162         }
03163         
03164         char result[10] ;
03165         size_t strIndx = 0;
03166         
03167         if (chargeFields) {
03168             while ((strIndx < (strlen(line)-9)) && (strlen(line)-1 >=9 ) ) {
03169                 
03170                 strncpy(result, line+strIndx,9) ;
03171                 result[9] = '\0';
03172                 
03173                 Real localCharge = atof(result);
03174                 
03175                 // If we are reading charges from QM atoms, store them.
03176                 if (atmIndx < msg->numQMAtoms ) {
03177                     atmP[atmIndx].charge = localCharge;
03178                     resForce[atmIndx].charge = localCharge;
03179                 }
03180                 
03181                 // If we are reading charges from Dummy atoms,
03182                 // place them on the appropriate QM atoms.
03183                 if ( msg->numQMAtoms <= atmIndx &&
03184                     atmIndx < msg->numAllAtoms ) {
03185                     // We keep the calculated charge in the dummy atom
03186                     // so that we can calculate electrostatic forces later on.
03187                     atmP[atmIndx].charge = localCharge;
03188                     
03189                     // The dummy atom points to the QM atom to which it is bound.
03190                     int qmInd = atmP[atmIndx].bountToIndx ;
03191                     resForce[qmInd].charge += localCharge;
03192                 }
03193                 
03194                 strIndx += 9;
03195                 atmIndx++;
03196                 
03197                 // If we found all charges for QM and dummy atoms, break the loop
03198                 // and stop reading the line.
03199                 if (atmIndx == msg->numAllAtoms) {
03200                     chargeFields = false;
03201                     break;
03202                 }
03203                 
03204             }
03205             
03206         }
03207         
03208         if (gradFields) {
03209             while ((strIndx < (strlen(line)-9)) && (strlen(line)-1 >=9 ) ) {
03210                 
03211                 strncpy(result, line+strIndx,9) ;
03212                 result[9] = '\0';
03213                 
03214                 gradient[gradCount] = atof(result);
03215                 if (gradCount == 2) {
03216                     
03217                     if (atmIndx < msg->numQMAtoms ) {
03218                         // Gradients in MOPAC are written in kcal/mol/A.
03219                         resForce[atmIndx].force.x = -1*gradient[0];
03220                         resForce[atmIndx].force.y = -1*gradient[1];
03221                         resForce[atmIndx].force.z = -1*gradient[2];
03222                     }
03223                     
03224                     // If we are reading forces applied on Dummy atoms,
03225                     // place them on the appropriate QM and MM atoms to conserve energy.
03226                     
03227                     // This implementation was based on the description in 
03228                     // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
03229                     if ( msg->numQMAtoms <= atmIndx &&
03230                     atmIndx < msg->numAllAtoms ) {
03231                         // The dummy atom points to the QM atom to which it is bound.
03232                         // The QM atom and the MM atom (in a QM-MM bond) point to each other.
03233                         int qmInd = atmP[atmIndx].bountToIndx ;
03234                         int mmInd = atmP[qmInd].bountToIndx ;
03235                         
03236                         Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
03237                         Real mmqmDist = dir.length() ;
03238                         
03239                         Real linkDist = Vector(atmP[atmIndx].position - 
03240                                         atmP[qmInd].position).length() ;
03241                         
03242                         Force mmForce(0), qmForce(0), 
03243                             linkForce(gradient[0], gradient[1], gradient[2]);
03244                         linkForce *= -1;
03245                         
03246                         Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
03247                         // Unit vectors
03248                         Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
03249                         Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
03250                         Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
03251                         Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
03252                         
03253                         qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv + 
03254                                     (xDelta)*base) )*xuv;
03255                         
03256                         qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv + 
03257                                     (yDelta)*base) )*yuv;
03258                         
03259                         qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv + 
03260                                     (zDelta)*base) )*zuv;
03261                         
03262                         
03263                         mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
03264                                     (xDelta)*base) )*xuv;
03265                         
03266                         mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
03267                                     (yDelta)*base) )*yuv;
03268                         
03269                         mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
03270                                     (zDelta)*base) )*zuv;
03271                         
03272                         resForce[qmInd].force += qmForce;
03273                         resForce[msg->numQMAtoms + mmInd].force += mmForce;
03274                     }
03275                     
03276                     gradCount = 0;
03277                     atmIndx++;
03278                 } else {
03279                     gradCount++;
03280                 }
03281                 
03282                 strIndx += 9;
03283                 
03284                 // If we found all gradients for QM atoms, break the loop
03285                 // and keep the next gradient line from being read, if there
03286                 // is one. Following gradients, if present, will refer to link
03287                 // atoms, and who cares about those?.
03288                 if (atmIndx == msg->numAllAtoms) {
03289                     gradFields = false;
03290                     break;
03291                 }
03292             }
03293         }
03294         
03295     }
03296     
03297     delete [] line;
03298     
03299     fclose(outputFile);
03300     
03301     // In case charges are not to be read form the QM software output,
03302     // we load the origianl atom charges.
03303     if (msg->qmAtmChrgMode == QMCHRGNONE) {
03304         int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
03305         const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
03306         Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
03307         
03308         atmIndx = 0 ;
03309         for (; atmIndx < msg->numQMAtoms; atmIndx++) {
03310             
03311             // Loops over all QM atoms (in all QM groups) comparing their global indices
03312             for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
03313                 
03314                 if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
03315                     
03316                     atmP[atmIndx].charge = qmAtmChrg[qmIter];
03317                     resForce[atmIndx].charge = qmAtmChrg[qmIter];
03318                     
03319                     break;
03320                 }
03321                 
03322             }
03323             
03324         }
03325         for (size_t j=msg->numQMAtoms; j<msg->numAllAtoms; ++j ) {
03326             atmP[j].charge = 0;
03327         }
03328     }
03329     
03330     // remove force file
03331 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
03332 //     iret = remove(outputFileName);
03333 //     if ( iret ) { NAMD_err(0); }
03334     
03335     if (! (chargesRead && gradsRead) ) {
03336         NAMD_die("Error reading QM forces file. Not all data could be read!");
03337     }
03338     
03339     DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
03340     
03341     atmP = msg->data ;
03342     pcP = msg->data + msg->numAllAtoms ;
03343     
03344     // The initial point charge index for the *force* message is the number of QM
03345     // atoms, since the dummy atoms have no representation in NAMD
03346     int pcIndx = msg->numQMAtoms;
03347     
03348     // ---- WARNING  ----
03349             // We need to apply the force felt by each QM atom due to the classical 
03350             // point charges sent to MOPAC.
03351             // MOPAC takes the MM electrostatic potential into cosideration ONLY
03352             // when performing the SCF calculation and when returning the energy
03353             // of the system, but it does not apply the potential to the final
03354             // gradient calculation, therefore, we calculate the resulting force
03355             // here.
03356     // Initialize force vector for electrostatic contribution
03357     std::vector<Force> qmElecForce ;
03358     for (size_t j=0; j<msg->numAllAtoms; ++j )
03359         qmElecForce.push_back( Force(0) );
03360     
03361     // We loop over point charges and distribute the forces applied over
03362     // virtual point charges to the MM1 and MM2 atoms (if any virtual PCs exists).
03363     for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
03364         
03365         // No force was applied to the QM region due to this charge, so it 
03366         // does not receive any back from the QM region. It must be an MM1 atom
03367         // from a QM-MM bond.
03368         if (pcP[i].type == QMPCTYPE_IGNORE)
03369             continue;
03370         
03371         Force totalForce(0);
03372         
03373         BigReal pntCharge = pcP[i].charge;
03374         
03375         Position posMM = pcP[i].position ;
03376         
03377         for (size_t j=0; j<msg->numAllAtoms; ++j ) {
03378             
03379             BigReal qmCharge = atmP[j].charge ;
03380             
03381             BigReal force = pntCharge*qmCharge*constants ;
03382             
03383             Position rVec = posMM - atmP[j].position ;
03384             
03385             force /= rVec.length2();
03386             
03387             // We accumulate the total force felt by a point charge
03388             // due to the QM charge distribution. This total force
03389             // will be applied to the point charge if it is a real one,
03390             // or will be distirbuted over MM1 and MM2 point charges, it 
03391             // this is a virtual point charge.
03392             totalForce += force*rVec.unit();
03393             
03394             // Accumulate force felt by a QM atom due to point charges.
03395             qmElecForce[j] += -1*force*rVec.unit();
03396         }
03397         
03398         if (pcP[i].type == QMPCTYPE_CLASSICAL) {
03399             // If we already ignored MM1 charges, we take all other 
03400             // non-virtual charges and apply forces directly to them.
03401             resForce[pcIndx].force += totalForce;
03402         }
03403         else {
03404             // If we are handling virtual PC, we distribute the force over
03405             // MM1 and MM2.
03406             
03407             // Virtual PC are bound to MM2.
03408             int mm2Indx = pcP[i].bountToIndx;
03409             // MM2 charges are bound to MM1.
03410             int mm1Indx = pcP[mm2Indx].bountToIndx;
03411             
03412             Real Cq = pcP[i].dist;
03413             
03414             Force mm1Force = (1-Cq)*totalForce ;
03415             Force mm2Force = Cq*totalForce ;
03416             
03417             resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
03418             resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
03419         }
03420         
03421     }
03422     
03423     // We now apply the accumulated electrostatic force contribution
03424     // to QM atoms.
03425     for (size_t j=0; j<msg->numAllAtoms; ++j ) {
03426         
03427         if (j < msg->numQMAtoms ) {
03428             resForce[j].force += qmElecForce[j];
03429             
03430         } else {
03431             // In case the QM atom is a Link atom, we redistribute 
03432             // its force as bove.
03433             
03434             // The dummy atom points to the QM atom to which it is bound.
03435             // The QM atom and the MM atom (in a QM-MM bond) point to each other.
03436             int qmInd = atmP[j].bountToIndx ;
03437             int mmInd = atmP[qmInd].bountToIndx ;
03438             
03439             Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
03440             Real mmqmDist = dir.length() ;
03441             
03442             Real linkDist = Vector(atmP[j].position - 
03443                             atmP[qmInd].position).length() ;
03444             
03445             Force linkForce = qmElecForce[j];
03446             
03447             Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
03448             // Unit vectors
03449             Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
03450             Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
03451             Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
03452             Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
03453             
03454             Force qmForce = (linkForce*((1 - linkDist/mmqmDist)*xuv + 
03455                         (xDelta)*base) )*xuv;
03456             
03457             qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv + 
03458                         (yDelta)*base) )*yuv;
03459             
03460             qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv + 
03461                         (zDelta)*base) )*zuv;
03462             
03463             
03464             Force mmForce = (linkForce*((linkDist/mmqmDist)*xuv -
03465                         (xDelta)*base) )*xuv;
03466             
03467             mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
03468                         (yDelta)*base) )*yuv;
03469             
03470             mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
03471                         (zDelta)*base) )*zuv;
03472             
03473             resForce[qmInd].force += qmForce;
03474             resForce[msg->numQMAtoms + mmInd].force += mmForce;
03475         }
03476     }
03477     
03478     // Adjusts forces from PME, canceling contributions from the QM and 
03479     // direct Coulomb forces calculated here.
03480     if (msg->PMEOn) {
03481         
03482         DebugM(1,"Correcting forces and energy for PME.\n");
03483         
03484         ewaldcof = msg->PMEEwaldCoefficient;
03485         BigReal TwoBySqrtPi = 1.12837916709551;
03486         pi_ewaldcof = TwoBySqrtPi * ewaldcof;
03487         
03488         for (size_t i=0; i < msg->numQMAtoms; i++) {
03489             
03490             BigReal p_i_charge = atmP[i].charge ;
03491             Position pos_i = atmP[i].position ;
03492             
03493             const BigReal kq_i = p_i_charge * constants;
03494             
03495             for (size_t j=i+1; j < msg->numQMAtoms; j++) {
03496                 
03497                 BigReal p_j_charge = atmP[j].charge ;
03498                 
03499                 Position pos_j = atmP[j].position ;
03500                 
03501                 BigReal r = Vector(pos_i - pos_j).length() ;
03502                 
03503                 BigReal tmp_a = r * ewaldcof;
03504                 BigReal tmp_b = erfc(tmp_a);
03505                 BigReal corr_energy = tmp_b;
03506                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
03507                 
03508 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
03509                 BigReal recip_energy = (1-tmp_b)/r;
03510                 
03511 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
03512                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
03513                 
03514                 // Final force and energy correction for this pair of atoms.
03515                 BigReal energy = kq_i * p_j_charge * recip_energy ;
03516                 
03517                 Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
03518                 
03519                 // The force is *subtracted* from the total force acting on
03520                 // both atoms. The sign on fixForce corrects the orientation
03521                 // of the subtracted force.
03522 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
03523 //                     << std::endl);
03524 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
03525 //                     << std::endl);
03526 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
03527 //                 DebugM(4,"Energy correction: " << energy << "\n");
03528                 resForce[i].force -= fixForce ;
03529                 resForce[j].force -= -1*fixForce ;
03530                 
03531                 // The energy is *subtracted* from the total energy calculated here.
03532                 resMsg->energyCorr -= energy;
03533             }
03534         }
03535         
03536 //         DebugM(4,"Corrected energy QM-QM interactions: " << resMsg->energyCorr << "\n");
03537         
03538         pcIndx = msg->numQMAtoms;
03539         // We only loop over point charges from real atoms, ignoring the ones 
03540         // created to handle QM-MM bonds.
03541         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
03542             
03543             BigReal p_i_charge = pcPpme[i].charge;
03544             Position pos_i = pcPpme[i].position ;
03545             
03546             const BigReal kq_i = p_i_charge * constants;
03547             
03548             Force fixForce = 0;
03549             
03550             for (size_t j=0; j<msg->numQMAtoms; ++j ) {
03551                 
03552                 BigReal p_j_charge = atmP[j].charge ;
03553                 
03554                 Position pos_j = atmP[j].position ;
03555                 
03556                 BigReal r = Vector(pos_i - pos_j).length() ;
03557                 
03558                 BigReal tmp_a = r * ewaldcof;
03559                 BigReal tmp_b = erfc(tmp_a);
03560                 BigReal corr_energy = tmp_b;
03561                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
03562                 
03563 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
03564                 BigReal recip_energy = (1-tmp_b)/r;
03565                 
03566 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
03567                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
03568                 
03569                 // Final force and energy correction for this pair of atoms.
03570                 BigReal energy = kq_i * p_j_charge * recip_energy ;
03571                 
03572                 fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
03573                 
03574                 // The energy is *subtracted* from the total energy calculated here.
03575                 resMsg->energyCorr -= energy;
03576             }
03577             
03578             // The force is *subtracted* from the total force acting on
03579             // the point charge..
03580 //             DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
03581 //             << std::endl);
03582 //             DebugM(4,"Force correction: " << fixForce << std::endl);
03583             resForce[pcIndx].force -= kq_i*fixForce ;
03584         }
03585         
03586     }
03587     
03588     // Calculates the virial contribution form the forces on QM atoms and 
03589     // point charges calculated here.
03590     for (size_t i=0; i < msg->numQMAtoms; i++) {
03591         // virial used by NAMD is -'ve of normal convention, so reverse it!
03592         // virial[i][j] in file should be sum of -1 * f_i * r_j
03593         for ( int k=0; k<3; ++k )
03594             for ( int l=0; l<3; ++l )
03595                 resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
03596     }
03597     
03598     pcIndx = msg->numQMAtoms; // Index in the real PC force array.
03599     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
03600         // virial used by NAMD is -'ve of normal convention, so reverse it!
03601         // virial[i][j] in file should be sum of -1 * f_i * r_j
03602         for ( int k=0; k<3; ++k )
03603             for ( int l=0; l<3; ++l )
03604                 resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
03605     }
03606     
03607     
03608     
03609     // Send message to rank zero with results.
03610     QMProxy[0].recvQMRes(resMsg);
03611     
03612     if (msg->switching && msg->PMEOn)
03613         delete [] pcPpme;
03614     
03615     delete msg;
03616     return ;
03617 }
03618 
03619 void ComputeQMMgr::calcORCA(QMGrpCalcMsg *msg)
03620 {
03621     
03622     FILE *inputFile,*outputFile,*chrgFile;
03623     int iret;
03624     
03625     const size_t lineLen = 256;
03626     char *line = new char[lineLen];
03627     
03628     std::string qmCommand, inputFileName, outputFileName, pntChrgFileName, pcGradFilename;
03629     std::string tmpRedirectFileName, pcGradFileName;
03630     
03631     // For coulomb calculation
03632     BigReal constants = msg->constants;
03633     
03634     double gradient[3];
03635     
03636     DebugM(4,"Running ORCA on PE " << CkMyPe() << std::endl);
03637     
03638     QMAtomData *pcP = msg->data + msg->numAllAtoms ;
03639     
03640     // We create a pointer for PME correction, which may point to
03641     // a copy of the original point charge array, with unchanged charges,
03642     // or points to the original array in case no switching or charge
03643     // scheme is being used.
03644     QMAtomData *pcPpme = NULL;
03645     if (msg->switching) {
03646         
03647         if (msg->PMEOn)
03648             pcPpme = new QMAtomData[msg->numRealPntChrgs];
03649         
03650         pntChrgSwitching(msg, pcPpme) ;
03651     } else {
03652         pcPpme = pcP;
03653     }
03654     
03655     // For each QM group, create a subdirectory where files will be palced.
03656     std::string baseDir(msg->baseDir);
03657     baseDir.append("/") ;
03658     std::ostringstream itosConv ;
03659     itosConv << msg->peIter ;
03660     baseDir += itosConv.str() ;
03661     
03662     struct stat info;
03663     
03664     if (stat(msg->baseDir, &info) != 0 ) {
03665         CkPrintf( "Node %d cannot access directory %s\n",
03666                   CkMyPe(), baseDir.c_str() );
03667         NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
03668     }
03669     else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
03670         DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
03671         int retVal = mkdir(baseDir.c_str(), S_IRWXU);
03672     }
03673     
03674     #ifdef DEBUG_QM
03675     Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
03676     #endif
03677     
03678     inputFileName.clear();
03679     inputFileName.append(baseDir.c_str()) ;
03680     inputFileName.append("/qmmm_") ;
03681     inputFileName += itosConv.str() ;
03682     inputFileName.append(".input") ;
03683     
03684     // Opens file for coordinate and parameter input
03685     inputFile = fopen(inputFileName.c_str(),"w");
03686     if ( ! inputFile ) {
03687         iout << iERROR << "Could not open input file for writing: " 
03688         << inputFileName << "\n" << endi ;
03689         NAMD_err("Error writing QM input file.");
03690     }
03691     
03692     // Builds the command that will be executed
03693     qmCommand.clear();
03694     qmCommand.append("cd ");
03695     qmCommand.append(baseDir);
03696     qmCommand.append(" ; ");
03697     qmCommand.append(msg->execPath) ;
03698     qmCommand.append(" ") ;
03699     qmCommand.append(inputFileName) ;
03700     
03701     // Build the redirect file name we need for the screen output.
03702     // That's the only place where we can get partial charges for QM atoms.
03703     tmpRedirectFileName = inputFileName ;
03704     tmpRedirectFileName.append(".TmpOut") ;
03705     
03706     qmCommand.append(" > ") ;
03707     qmCommand.append(tmpRedirectFileName) ;
03708     
03709     // Builds the file name where orca will place the gradient
03710     // This will be relative to the input file
03711     outputFileName = inputFileName ;
03712     outputFileName.append(".engrad") ;
03713     
03714     // Builds the file name where orca will place gradients acting on 
03715     // point charges
03716     pcGradFilename = inputFileName ;
03717     pcGradFilename.append(".pcgrad") ;
03718     
03719     if (msg->numAllPntChrgs) {
03720         // Builds the file name where we will write the point charges.
03721         pntChrgFileName = inputFileName ;
03722         pntChrgFileName.append(".pntchrg") ;
03723         
03724         pcGradFileName = inputFileName;
03725         pcGradFileName.append(".pcgrad");
03726         
03727         chrgFile = fopen(pntChrgFileName.c_str(),"w");
03728         if ( ! chrgFile ) {
03729             iout << iERROR << "Could not open charge file for writing: " 
03730             << pntChrgFileName << "\n" << endi ;
03731             NAMD_err("Error writing point charge file.");
03732         }
03733         
03734         int numPntChrgs = 0;
03735         for (int i=0; i<msg->numAllPntChrgs; i++ ) {
03736             if (pcP[i].type != QMPCTYPE_IGNORE)
03737                 numPntChrgs++;
03738         }
03739         
03740         iret = fprintf(chrgFile,"%d\n", numPntChrgs);
03741         if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
03742     }
03743     
03744     // writes configuration lines to the input file.
03745     std::stringstream ss(msg->configLines) ;
03746     std::string confLineStr;
03747     while (std::getline(ss, confLineStr) ) {
03748         confLineStr.append("\n");
03749         iret = fprintf(inputFile,confLineStr.c_str());
03750         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03751     }
03752     
03753     if (msg->numAllPntChrgs) {
03754         iret = fprintf(inputFile,"%%pointcharges \"%s\"\n", pntChrgFileName.c_str());
03755         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03756     }
03757     
03758     iret = fprintf(inputFile,"\n\n%%coords\n  CTyp xyz\n");
03759     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03760     
03761     iret = fprintf(inputFile,"  Charge %f\n",msg->charge);
03762     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03763     
03764     iret = fprintf(inputFile,"  Mult %f\n",msg->multiplicity);
03765     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03766     
03767     iret = fprintf(inputFile,"  Units Angs\n  coords\n\n");
03768     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03769     
03770     DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " << 
03771     inputFileName.c_str() << " and " << msg->numAllPntChrgs << 
03772     " point charges in file " << pntChrgFileName.c_str() << "\n");
03773     
03774     // write QM and dummy atom coordinates to input file.
03775     QMAtomData *atmP = msg->data ;
03776     for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
03777         
03778         double x = atmP->position.x;
03779         double y = atmP->position.y;
03780         double z = atmP->position.z;
03781         
03782         iret = fprintf(inputFile,"  %s %f %f %f\n",
03783                        atmP->element,x,y,z);
03784         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03785         
03786     }
03787     
03788     iret = fprintf(inputFile,"  end\nend\n");
03789     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03790     
03791     if (msg->numAllPntChrgs) {
03792         // Write point charges to file.
03793         pcP = msg->data + msg->numAllAtoms ;
03794         for ( size_t j=0; j < msg->numAllPntChrgs; j++, ++pcP) {
03795             
03796             if (pcP->type == QMPCTYPE_IGNORE)
03797                 continue;
03798             
03799             double charge = pcP->charge;
03800             
03801             double x = pcP->position.x;
03802             double y = pcP->position.y;
03803             double z = pcP->position.z;
03804             
03805             iret = fprintf(chrgFile,"%f %f %f %f\n",
03806                            charge,x,y,z);
03807             if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
03808         }
03809         
03810         fclose(chrgFile);
03811     }
03812     
03813     DebugM(4,"Closing input and charge file\n");
03814     fclose(inputFile);
03815     
03816     if (msg->prepProcOn) {
03817         
03818         std::string prepProc(msg->prepProc) ;
03819         prepProc.append(" ") ;
03820         prepProc.append(inputFileName) ;
03821         iret = system(prepProc.c_str());
03822         if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
03823         if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
03824     }
03825     
03826         // runs QM command
03827     DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
03828     iret = system(qmCommand.c_str());
03829     
03830     if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
03831     if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
03832 
03833     if (msg->secProcOn) {
03834         
03835         std::string secProc(msg->secProc) ;
03836         secProc.append(" ") ;
03837         secProc.append(inputFileName) ;
03838         itosConv.str("");
03839         itosConv.clear() ;
03840         itosConv << msg->timestep ;
03841         secProc.append(" ") ;
03842         secProc += itosConv.str() ;
03843         
03844         iret = system(secProc.c_str());
03845         if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
03846         if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
03847     }
03848 
03849     // remove coordinate file
03850 //     iret = remove(inputFileName);
03851 //     if ( iret ) { NAMD_err(0); }
03852 
03853     // remove coordinate file
03854 //     iret = remove(pntChrgFileName);
03855 //     if ( iret ) { NAMD_err(0); }
03856     
03857     // opens output file
03858     DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
03859     outputFile = fopen(outputFileName.c_str(),"r");
03860     if ( ! outputFile ) {
03861         iout << iERROR << "Could not find QM output file!\n" << endi;
03862         NAMD_err(0); 
03863     }
03864 
03865     // Resets the pointers.
03866     atmP = msg->data ;
03867     pcP = msg->data + msg->numAllAtoms ;
03868     
03869     // Allocates the return message.
03870     QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
03871     resMsg->grpIndx = msg->grpIndx;
03872     resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
03873     resMsg->energyOrig = 0;
03874     resMsg->energyCorr = 0;
03875     for ( int k=0; k<3; ++k )
03876         for ( int l=0; l<3; ++l )
03877             resMsg->virial[k][l] = 0;
03878     QMForce *resForce = resMsg->force ;
03879     
03880     // We split the data into two pointers so we can "skip" the dummy atoms
03881     // which have no representation in NAMD.
03882     for (int i=0; i<resMsg->numForces; i++) {
03883         resForce[i].force = 0;
03884         resForce[i].charge = 0 ;
03885         if (i < msg->numQMAtoms) {
03886             // We use the replace field to indicate QM atoms,
03887             // which will have charge information.
03888             resForce[i].replace = 1 ;
03889             resForce[i].id = atmP->id;
03890             atmP++;
03891         }
03892         else {
03893             // We use the replace field to indicate QM atoms,
03894             // which will have charge information.
03895             resForce[i].replace = 0 ;
03896             resForce[i].id = pcP->id;
03897             pcP++;
03898         }
03899     }
03900     
03901     // Resets the pointers.
03902     atmP = msg->data ;
03903     pcP = msg->data + msg->numAllAtoms ;
03904     
03905     size_t atmIndx = 0, gradCount = 0;
03906     int numAtomsInFile = 0;
03907     
03908     // Reads the data form the output file created by the QM software.
03909     // Gradients over the QM atoms, and Charges for QM atoms will be read.
03910     
03911     // skip lines before number of atoms
03912     for (int i = 0; i < 3; i++) {
03913         fgets(line, lineLen, outputFile); 
03914     }
03915     
03916     iret = fscanf(outputFile,"%d\n", &numAtomsInFile);
03917     if ( iret != 1 ) {
03918         NAMD_die("Error reading QM forces file.");
03919     }
03920     
03921     #ifdef DEBUG_QM
03922     if(numAtomsInFile != msg->numAllAtoms) {
03923         NAMD_die("Error reading QM forces file. Number of atoms in QM output\
03924         does not match the expected.");
03925     }
03926     #endif
03927     
03928     // skip lines before energy
03929     for (int i = 0; i < 3; i++) {
03930         fgets(line, lineLen, outputFile); 
03931     }
03932     
03933     iret = fscanf(outputFile,"%lf\n", &resMsg->energyOrig);
03934     if ( iret != 1 ) {
03935         NAMD_die("Error reading QM forces file.");
03936     }
03937 //     iout << "Energy in step (Hartree): " << resMsg->energyOrig << "\n" <<  endi;
03938     // All energies are given in Eh (Hartree)
03939     // NAMD needs energies in kcal/mol
03940     // The conversion factor is 627.509469
03941     resMsg->energyOrig *= 627.509469;
03942 //     iout << "Energy in step (Kcal/mol): " << resMsg->energyOrig << "\n" <<  endi;
03943     
03944     resMsg->energyCorr = resMsg->energyOrig;
03945     
03946     // skip lines before gradient
03947     for (int i = 0; i < 3; i++) {
03948         fgets(line, lineLen, outputFile) ;
03949     }
03950     
03951     // Break the loop when we find all gradients for QM atoms, 
03952     // and keep the next gradient lines from being read, if there
03953     // are more. Following gradients, if present, will refer to link
03954     // atoms.
03955     atmIndx = 0;
03956     gradCount = 0;
03957     for ( size_t i=0; i<msg->numAllAtoms*3; ++i ) {
03958         
03959         iret = fscanf(outputFile,"%lf\n", &gradient[gradCount]);
03960         if ( iret != 1 ) { NAMD_die("Error reading QM forces file."); }
03961         
03962         if (gradCount == 2){
03963             
03964             // All gradients are given in Eh/a0 (Hartree over Bohr radius)
03965             // NAMD needs forces in kcal/mol/angstrons
03966             // The conversion factor is 627.509469/0.529177 = 1185.82151
03967             if (atmIndx < msg->numQMAtoms ) {
03968                 resForce[atmIndx].force.x = -1*gradient[0]*1185.82151;
03969                 resForce[atmIndx].force.y = -1*gradient[1]*1185.82151;
03970                 resForce[atmIndx].force.z = -1*gradient[2]*1185.82151;
03971             }
03972             
03973             // If we are reading forces applied on Dummy atoms,
03974             // place them on the appropriate QM and MM atoms to conserve energy.
03975             
03976             // This implementation was based on the description in 
03977             // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
03978             if ( msg->numQMAtoms <= atmIndx &&
03979             atmIndx < msg->numAllAtoms ) {
03980                 // The dummy atom points to the QM atom to which it is bound.
03981                 // The QM atom and the MM atom (in a QM-MM bond) point to each other.
03982                 int qmInd = atmP[atmIndx].bountToIndx ;
03983                 int mmInd = atmP[qmInd].bountToIndx ;
03984                 
03985                 Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
03986                 Real mmqmDist = dir.length() ;
03987                 
03988                 Real linkDist = Vector(atmP[atmIndx].position - atmP[qmInd].position).length() ;
03989                 
03990                 Force mmForce(0), qmForce(0), linkForce(gradient[0], gradient[1], gradient[2]);
03991                 linkForce *= -1*1185.82151;
03992                 
03993                 Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
03994                 // Unit vectors
03995                 Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
03996                 Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
03997                 Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
03998                 Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
03999                 
04000                 qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv + 
04001                             (xDelta)*base) )*xuv;
04002                 
04003                 qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv + 
04004                             (yDelta)*base) )*yuv;
04005                 
04006                 qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv + 
04007                             (zDelta)*base) )*zuv;
04008                 
04009                 
04010                 mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
04011                             (xDelta)*base) )*xuv;
04012                 
04013                 mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
04014                             (yDelta)*base) )*yuv;
04015                 
04016                 mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
04017                             (zDelta)*base) )*zuv;
04018                 
04019                 resForce[qmInd].force += qmForce;
04020                 resForce[msg->numQMAtoms + mmInd].force += mmForce;
04021             }
04022             
04023             gradCount = 0;
04024             atmIndx++ ;
04025         } else
04026             gradCount++ ;
04027         
04028     }
04029     
04030     fclose(outputFile);
04031     
04032     // In case charges are not to be read form the QM software output,
04033     // we load the origianl atom charges.
04034     if (msg->qmAtmChrgMode == QMCHRGNONE) {
04035         int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
04036         const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
04037         Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
04038         
04039         atmIndx = 0 ;
04040         for (; atmIndx < msg->numQMAtoms; atmIndx++) {
04041             
04042             // Loops over all QM atoms (in all QM groups) comparing their global indices
04043             for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
04044                 
04045                 if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
04046                     
04047                     atmP[atmIndx].charge = qmAtmChrg[qmIter];
04048                     resForce[atmIndx].charge = qmAtmChrg[qmIter];
04049                     
04050                     break;
04051                 }
04052                 
04053             }
04054             
04055         }
04056     }
04057     else {
04058         // opens redirected output file
04059         outputFile = fopen(tmpRedirectFileName.c_str(),"r");
04060         if ( ! outputFile ) {
04061             iout << iERROR << "Could not find Redirect output file:"
04062             << tmpRedirectFileName << std::endl << endi;
04063             NAMD_err(0); 
04064         }
04065         
04066         DebugM(4,"Opened tmeporary output for charge reading: " << tmpRedirectFileName.c_str() << "\n");
04067         
04068         int lineState = 0;
04069         char result[20] ;
04070         while ( fgets(line, lineLen, outputFile) != NULL){
04071             
04072             // We loop over the lines of the output file untill we find
04073             // the line that initiates the charges lines. Then
04074             // we read all lines and expect to get one charge
04075             // per line, untill we find a line that sums all charges.
04076             
04077             switch (msg->qmAtmChrgMode) {
04078                 
04079                 case QMCHRGMULLIKEN:
04080                 {
04081                 
04082                 if ( strstr(line,"MULLIKEN ATOMIC CHARGES") != NULL ) {
04083                     lineState = 1;
04084                     atmIndx = 0;
04085                     
04086                     // Now we expect the following line to have a series of dashes
04087                     // and the folowing lines to have atom charges (among other info)
04088                     continue;
04089                 }
04090                 
04091                 if ( strstr(line,"Sum of atomic charges") != NULL ) {
04092                     
04093                     // Now that all charges have been read, grabs the total charge
04094                     // just for fun... and checking purposes.
04095                     #ifdef DEBUG_QM
04096                     strncpy(result, line + 31,12) ;
04097                     result[12] = '\0';
04098                     
04099                     DebugM(4,"Total charge of QM region calculated by ORCA is: " 
04100                     << atof(result) << std::endl )
04101                     #endif
04102                     
04103                     // Nothing more to read from the output file
04104                     break;
04105                 }
04106                 
04107                 // Line state 1 means we have to skip ONLY one line.
04108                 if (lineState == 1) {
04109                     lineState = 2;
04110                     continue;
04111                 }
04112                 
04113                 // Line state 2 means we have to read the line and grab the info.
04114                 if (lineState == 2) {
04115                     
04116                     strncpy(result, line+8,12) ;
04117                     result[12] = '\0';
04118                     
04119                     Real localCharge = atof(result);
04120                     
04121                     // If we are reading charges from QM atoms, store them.
04122                     if (atmIndx < msg->numQMAtoms ) {
04123                         atmP[atmIndx].charge = localCharge;
04124                         resForce[atmIndx].charge = localCharge;
04125                     }
04126                     
04127                     // If we are reading charges from Dummy atoms,
04128                     // place the on the appropriate QM atom.
04129                      if ( msg->numQMAtoms <= atmIndx &&
04130                         atmIndx < msg->numAllAtoms ) {
04131                         int qmInd = atmP[atmIndx].bountToIndx ;
04132                         atmP[qmInd].charge += localCharge;
04133                         resForce[qmInd].charge += localCharge;
04134                     }
04135                     
04136                     atmIndx++ ;
04137                     
04138                     // If we found all charges for QM atoms, change the lineState
04139                     // untill we reach the "total charge" line.
04140                     if (atmIndx == msg->numAllAtoms ) {
04141                         lineState = 0;
04142                     }
04143                     
04144                     continue;
04145                 }
04146                 
04147                 } break ;
04148                 
04149                 case QMCHRGCHELPG :
04150                 {
04151                 
04152                 if ( strstr(line,"CHELPG Charges") != NULL ) {
04153                     lineState = 1;
04154                     atmIndx = 0;
04155                     
04156                     // Now we expect the following line to have a series of dashes
04157                     // and the folowing lines to have atom charges (among other info)
04158                     continue;
04159                 }
04160                 
04161                 if ( strstr(line,"Total charge") != NULL ) {
04162                     
04163                     // Now that all charges have been read, grabs the total charge
04164                     // just for fun... and checking purposes.
04165                     #ifdef DEBUG_QM
04166                     strncpy(result, line + 14,13) ;
04167                     result[13] = '\0';
04168                     
04169                     DebugM(4,"Total charge of QM region calculated by ORCA is: " 
04170                     << atof(result) << std::endl )
04171                     #endif
04172                     
04173                     // Nothing more to read from the output file
04174                     break;
04175                 }
04176                 
04177                 // Line state 1 means we have to skip ONLY one line.
04178                 if (lineState == 1) {
04179                     lineState = 2;
04180                     continue;
04181                 }
04182                 
04183                 // Line state 2 means we have to read the line and grab the info.
04184                 if (lineState == 2) {
04185                     
04186                     strncpy(result, line+12,15) ;
04187                     result[15] = '\0';
04188                     
04189                     Real localCharge = atof(result);
04190                     
04191                     // If we are reading charges from QM atoms, store them.
04192                     if (atmIndx < msg->numQMAtoms ) {
04193                         atmP[atmIndx].charge = localCharge;
04194                         resForce[atmIndx].charge = localCharge;
04195                     }
04196                     
04197                     // If we are reading charges from Dummy atoms,
04198                     // place the on the appropriate QM atom.
04199                      if ( msg->numQMAtoms <= atmIndx &&
04200                         atmIndx < msg->numAllAtoms ) {
04201                         int qmInd = atmP[atmIndx].bountToIndx ;
04202                         atmP[qmInd].charge += localCharge;
04203                         resForce[qmInd].charge += localCharge;
04204                     }
04205                     
04206                     atmIndx++ ;
04207                     
04208                     // If we found all charges for QM atoms, we ignore the following line
04209                     // untill we reach the "total charge" line.
04210                     if (atmIndx == msg->numAllAtoms ) {
04211                         lineState = 1;
04212                     }
04213                     
04214                     continue;
04215                 }
04216                 
04217                 } break;
04218             }
04219         }
04220         
04221         fclose(outputFile);
04222     }
04223     
04224     // remove force file
04225 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
04226 //     iret = remove(outputFileName);
04227 //     if ( iret ) { NAMD_err(0); }
04228     
04229     
04230     DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
04231     
04232     atmP = msg->data ;
04233     pcP = msg->data + msg->numAllAtoms ;
04234     
04235     // The initial point charge index for the force message is the number of QM
04236     // atoms, since the dummy atoms have no representation in NAMD
04237     int pcIndx = msg->numQMAtoms;
04238     
04239     // If there are no point charges, orca will not create a "pcgrad" file.
04240     if (msg->numAllPntChrgs > 0) {
04241     
04242         // opens redirected output file
04243         outputFile = fopen(pcGradFileName.c_str(),"r");
04244         if ( ! outputFile ) {
04245             iout << iERROR << "Could not find PC gradient output file:"
04246             << pcGradFileName << std::endl << endi;
04247             NAMD_err(0); 
04248         }
04249         
04250         DebugM(4,"Opened pc output for gradient reading: " << pcGradFileName.c_str() << "\n");
04251         
04252         int pcgradNumPC = 0, readPC = 0;
04253         iret = fscanf(outputFile,"%d\n", &pcgradNumPC );
04254         if ( iret != 1 ) {
04255             NAMD_die("Error reading PC forces file.");
04256         }
04257         DebugM(4,"Found in pc gradient output: " << pcgradNumPC << " point charge grads.\n");
04258         
04259         // We loop over point charges, reading the total electrostatic force 
04260         // applied on them by the QM region.
04261         // We redistribute the forces applied over virtual point
04262         // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
04263         for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
04264             
04265             Force totalForce(0);
04266             
04267             // No force was applied to the QM region due to this charge, since it
04268             // was skipped when writing down point charges to the QM software, so it
04269             // does not receive any back from the QM region. It must be an MM1 atom
04270             // from a QM-MM bond.
04271             if (pcP[i].type == QMPCTYPE_IGNORE)
04272                 continue;
04273             
04274             fgets(line, lineLen, outputFile);
04275             
04276             iret = sscanf(line,"%lf %lf %lf\n", &totalForce[0], &totalForce[1], &totalForce[2] );
04277             if ( iret != 3 ) {
04278                 NAMD_die("Error reading PC forces file.");
04279             }
04280             // All gradients are given in Eh/a0 (Hartree over Bohr radius)
04281             // NAMD needs forces in kcal/mol/angstrons
04282             // The conversion factor is 627.509469/0.529177 = 1185.82151
04283             totalForce *= -1185.82151;
04284             
04285             if (pcP[i].type == QMPCTYPE_CLASSICAL) {
04286                 // If we already ignored MM1 charges, we take all other 
04287                 // non-virtual charges and apply forces directly to them.
04288                 resForce[pcIndx].force += totalForce;
04289             }
04290             else {
04291                 // If we are handling virtual PC, we distribute the force over
04292                 // MM1 and MM2.
04293                 
04294                 Force mm1Force(0), mm2Force(0);
04295                 
04296                 // Virtual PC are bound to MM2.
04297                 int mm2Indx = pcP[i].bountToIndx;
04298                 // MM2 charges are bound to MM1.
04299                 int mm1Indx = pcP[mm2Indx].bountToIndx;
04300                 
04301                 Real Cq = pcP[i].dist;
04302                 
04303                 mm1Force = (1-Cq)*totalForce ;
04304                 mm2Force = Cq*totalForce ;
04305                 
04306                 resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
04307                 resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
04308             }
04309             
04310         }
04311         
04312         fclose(outputFile);
04313     }
04314     
04315     delete [] line;
04316     
04317     // Adjusts forces from PME, canceling contributions from the QM and 
04318     // direct Coulomb forces calculated here.
04319     if (msg->PMEOn) {
04320         
04321         DebugM(1,"Correcting forces and energy for PME.\n");
04322         
04323         ewaldcof = msg->PMEEwaldCoefficient;
04324         BigReal TwoBySqrtPi = 1.12837916709551;
04325         pi_ewaldcof = TwoBySqrtPi * ewaldcof;
04326         
04327         for (size_t i=0; i < msg->numQMAtoms; i++) {
04328             
04329             BigReal p_i_charge = atmP[i].charge ;
04330             Position pos_i = atmP[i].position ;
04331             
04332             for (size_t j=i+1; j < msg->numQMAtoms; j++) {
04333                 
04334                 BigReal p_j_charge = atmP[j].charge ;
04335                 
04336                 Position pos_j = atmP[j].position ;
04337                 
04338                 BigReal r = Vector(pos_i - pos_j).length() ;
04339                 
04340                 BigReal tmp_a = r * ewaldcof;
04341                 BigReal tmp_b = erfc(tmp_a);
04342                 BigReal corr_energy = tmp_b;
04343                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
04344                 
04345 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
04346                 BigReal recip_energy = (1-tmp_b)/r;
04347                 
04348 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
04349                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
04350                 
04351                 const BigReal kq_i = p_i_charge * constants;
04352                 
04353                 // Final force and energy correction for this pair of atoms.
04354                 BigReal energy = kq_i * p_j_charge * recip_energy ;
04355                 
04356                 Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
04357                 
04358                 // The force is *subtracted* from the total force acting on
04359                 // both atoms. The sign on fixForce corrects the orientation
04360                 // of the subtracted force.
04361 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
04362 //                     << std::endl);
04363 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
04364 //                     << std::endl);
04365 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
04366                 resForce[i].force -= fixForce ;
04367                 resForce[j].force -= -1*fixForce;
04368                 
04369                 // The energy is *subtracted* from the total energy calculated here.
04370                 resMsg->energyCorr -= energy;
04371             }
04372         }
04373         
04374         pcIndx = msg->numQMAtoms;
04375         // We only loop over point charges from real atoms, ignoring the ones 
04376         // created to handle QM-MM bonds.
04377         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
04378             
04379             BigReal p_i_charge = pcPpme[i].charge;
04380             Position pos_i = pcPpme[i].position ;
04381             
04382             const BigReal kq_i = p_i_charge * constants;
04383             
04384             Force fixForce = 0;
04385             
04386             for (size_t j=0; j<msg->numQMAtoms; ++j ) {
04387                 
04388                 BigReal p_j_charge = atmP[j].charge ;
04389                 
04390                 Position pos_j = atmP[j].position ;
04391                 
04392                 BigReal r = Vector(pos_i - pos_j).length() ;
04393                 
04394                 BigReal tmp_a = r * ewaldcof;
04395                 BigReal tmp_b = erfc(tmp_a);
04396                 BigReal corr_energy = tmp_b;
04397                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
04398                 
04399 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
04400                 BigReal recip_energy = (1-tmp_b)/r;
04401                 
04402 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
04403                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
04404                 
04405                 // Final force and energy correction for this pair of atoms.
04406                 BigReal energy = kq_i * p_j_charge * recip_energy ;
04407                 
04408                 fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
04409                 
04410                 // The energy is *subtracted* from the total energy calculated here.
04411                 resMsg->energyCorr -= energy;
04412                 
04413             }
04414             
04415             // The force is *subtracted* from the total force acting on
04416                 // the point charge.
04417 //                 DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
04418 //                     << std::endl);
04419 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
04420             resForce[pcIndx].force -= kq_i*fixForce ;
04421         }
04422         
04423     }
04424     
04425     DebugM(1,"Determining virial...\n");
04426     
04427     // Calculates the virial contribution form the forces on QM atoms and 
04428     // point charges calculated here.
04429     for (size_t i=0; i < msg->numQMAtoms; i++) {
04430         // virial used by NAMD is -'ve of normal convention, so reverse it!
04431         // virial[i][j] in file should be sum of -1 * f_i * r_j
04432         for ( int k=0; k<3; ++k )
04433             for ( int l=0; l<3; ++l )
04434                 resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
04435     }
04436     
04437     pcIndx = msg->numQMAtoms; // Index in the real PC force array.
04438     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
04439         // virial used by NAMD is -'ve of normal convention, so reverse it!
04440         // virial[i][j] in file should be sum of -1 * f_i * r_j
04441         for ( int k=0; k<3; ++k )
04442             for ( int l=0; l<3; ++l )
04443                 resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
04444     }
04445     
04446     DebugM(1,"End of QM processing. Sending result message.\n");
04447     
04448     // Send message to rank zero with results.
04449     QMProxy[0].recvQMRes(resMsg);
04450     
04451     if (msg->switching && msg->PMEOn)
04452         delete [] pcPpme;
04453     
04454     delete msg;
04455     return ;
04456 }
04457 
04458 void ComputeQMMgr::calcUSR(QMGrpCalcMsg *msg) {
04459     
04460     FILE *inputFile,*outputFile;
04461     int iret;
04462     
04463     std::string qmCommand, inputFileName, outputFileName;
04464     
04465     // For coulomb calculation
04466     BigReal constants = msg->constants;
04467     
04468     DebugM(4,"Running USER DEFINED SOFTWARE on PE " << CkMyPe() << std::endl);
04469     
04470     QMAtomData *pcP = msg->data + msg->numAllAtoms ;
04471     
04472     // We create a pointer for PME correction, which may point to
04473     // a copy of the original point charge array, with unchanged charges,
04474     // or points to the original array in case no switching or charge
04475     // scheme is being used.
04476     QMAtomData *pcPpme = NULL;
04477     if (msg->switching) {
04478         
04479         if (msg->PMEOn)
04480             pcPpme = new QMAtomData[msg->numRealPntChrgs];
04481         
04482         pntChrgSwitching(msg, pcPpme) ;
04483     } else {
04484         pcPpme = pcP;
04485     }
04486     
04487     // For each QM group, create a subdirectory where files will be palced.
04488     std::string baseDir(msg->baseDir);
04489     baseDir.append("/") ;
04490     std::ostringstream itosConv ;
04491     itosConv << msg->peIter ;
04492     baseDir += itosConv.str() ;
04493     
04494     struct stat info;
04495     
04496     if (stat(msg->baseDir, &info) != 0 ) {
04497         CkPrintf( "Node %d cannot access directory %s\n",
04498                   CkMyPe(), baseDir.c_str() );
04499         NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
04500     }
04501     else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
04502         DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
04503         int retVal = mkdir(baseDir.c_str(), S_IRWXU);
04504     }
04505     
04506     #ifdef DEBUG_QM
04507     Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
04508     #endif
04509     
04510     inputFileName.clear();
04511     inputFileName.append(baseDir.c_str()) ;
04512     inputFileName.append("/qmmm_") ;
04513     inputFileName += itosConv.str() ;
04514     inputFileName.append(".input") ;
04515     
04516     // Opens file for coordinate and parameter input
04517     inputFile = fopen(inputFileName.c_str(),"w");
04518     if ( ! inputFile ) {
04519         iout << iERROR << "Could not open input file for writing: " 
04520         << inputFileName << "\n" << endi ;
04521         NAMD_err("Error writing QM input file.");
04522     }
04523     
04524     // Builds the command that will be executed
04525     qmCommand.clear();
04526     qmCommand.append("cd ");
04527     qmCommand.append(baseDir);
04528     qmCommand.append(" ; ");
04529     qmCommand.append(msg->execPath) ;
04530     qmCommand.append(" ") ;
04531     qmCommand.append(inputFileName) ;
04532     
04533     // Builds the file name where orca will place the gradient
04534     // This will be relative to the input file
04535     outputFileName = inputFileName ;
04536     outputFileName.append(".result") ;
04537     
04538     int numPntChrgs = 0;
04539     for (int i=0; i<msg->numAllPntChrgs; i++ ) {
04540         if (pcP[i].type != QMPCTYPE_IGNORE)
04541             numPntChrgs++;
04542     }
04543     
04544     iret = fprintf(inputFile,"%d %d\n",msg->numAllAtoms, numPntChrgs);
04545     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
04546     
04547     DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " << 
04548         inputFileName.c_str() << " and " << msg->numAllPntChrgs << 
04549         " point charges." << std::endl);
04550     
04551     // write QM and dummy atom coordinates to input file.
04552     QMAtomData *atmP = msg->data ;
04553     for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
04554         
04555         double x = atmP->position.x;
04556         double y = atmP->position.y;
04557         double z = atmP->position.z;
04558         
04559         iret = fprintf(inputFile,"%f %f %f %s\n",
04560                        x,y,z,atmP->element);
04561         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
04562         
04563     }
04564     
04565     int numWritenPCs = 0;
04566     // Write point charges to file.
04567     pcP = msg->data + msg->numAllAtoms ;
04568     for ( size_t j=0; j < msg->numAllPntChrgs; j++, ++pcP) {
04569         
04570         if (pcP->type == QMPCTYPE_IGNORE)
04571                 continue;
04572         
04573         double charge = pcP->charge;
04574         
04575         double x = pcP->position.x;
04576         double y = pcP->position.y;
04577         double z = pcP->position.z;
04578         
04579         iret = fprintf(inputFile,"%f %f %f %f\n",
04580                        x,y,z,charge);
04581         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
04582         
04583         numWritenPCs++;
04584     }
04585     
04586     DebugM(4,"Closing input file\n");
04587     fclose(inputFile);
04588     
04589     if (msg->prepProcOn) {
04590         
04591         std::string prepProc(msg->prepProc) ;
04592         prepProc.append(" ") ;
04593         prepProc.append(inputFileName) ;
04594         iret = system(prepProc.c_str());
04595         if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
04596         if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
04597     }
04598     
04599         // runs QM command
04600     DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
04601     iret = system(qmCommand.c_str());
04602     
04603     if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
04604     if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
04605 
04606     if (msg->secProcOn) {
04607         
04608         std::string secProc(msg->secProc) ;
04609         secProc.append(" ") ;
04610         secProc.append(inputFileName) ;
04611         itosConv.str("");
04612         itosConv.clear() ;
04613         itosConv << msg->timestep ;
04614         secProc.append(" ") ;
04615         secProc += itosConv.str() ;
04616         
04617         iret = system(secProc.c_str());
04618         if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
04619         if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
04620     }
04621 
04622     // remove coordinate file
04623 //     iret = remove(inputFileName);
04624 //     if ( iret ) { NAMD_err(0); }
04625 
04626     // remove coordinate file
04627 //     iret = remove(pntChrgFileName);
04628 //     if ( iret ) { NAMD_err(0); }
04629     
04630     // opens output file
04631     DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
04632     outputFile = fopen(outputFileName.c_str(),"r");
04633     if ( ! outputFile ) {
04634         iout << iERROR << "Could not find QM output file!\n" << endi;
04635         NAMD_err(0); 
04636     }
04637 
04638     // Resets the pointers.
04639     atmP = msg->data ;
04640     pcP = msg->data + msg->numAllAtoms ;
04641     
04642     // Allocates the return message.
04643     QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
04644     resMsg->grpIndx = msg->grpIndx;
04645     resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
04646     resMsg->energyOrig = 0;
04647     resMsg->energyCorr = 0;
04648     for ( int k=0; k<3; ++k )
04649         for ( int l=0; l<3; ++l )
04650             resMsg->virial[k][l] = 0;
04651     QMForce *resForce = resMsg->force ;
04652     
04653     // We split the data into two pointers so we can "skip" the dummy atoms
04654     // which have no representation in NAMD.
04655     for (int i=0; i<resMsg->numForces; i++) {
04656         resForce[i].force = 0;
04657         resForce[i].charge = 0 ;
04658         if (i < msg->numQMAtoms) {
04659             // We use the replace field to indicate QM atoms,
04660             // which will have charge information.
04661             resForce[i].replace = 1 ;
04662             resForce[i].id = atmP->id;
04663             atmP++;
04664         }
04665         else {
04666             // We use the replace field to indicate QM atoms,
04667             // which will have charge information.
04668             resForce[i].replace = 0 ;
04669             resForce[i].id = pcP->id;
04670             pcP++;
04671         }
04672     }
04673     
04674     // Resets the pointers.
04675     atmP = msg->data ;
04676     pcP = msg->data + msg->numAllAtoms ;
04677     
04678     // Reads the data form the output file created by the QM software.
04679     // Gradients over the QM atoms, and Charges for QM atoms will be read.
04680     
04681     // Number of point charges for which we will receive forces.
04682     int usrPCnum = 0;
04683     const size_t lineLen = 256;
04684     char *line = new char[lineLen];
04685     
04686     fgets(line, lineLen, outputFile);
04687     
04688 //     iret = fscanf(outputFile,"%lf %d\n", &resMsg->energyOrig, &usrPCnum);
04689     iret = sscanf(line,"%lf %i\n", &resMsg->energyOrig, &usrPCnum);
04690     if ( iret < 1 ) {
04691         NAMD_die("Error reading energy from QM results file.");
04692     }
04693     
04694     resMsg->energyCorr = resMsg->energyOrig;
04695     
04696     if (iret == 2 && numWritenPCs != usrPCnum) {
04697         iout << iERROR << "Number of point charges does not match what was provided!\n" << endi ;
04698         NAMD_die("Error reading QM results file.");
04699     }
04700     
04701     size_t atmIndx;
04702     double localForce[3];
04703     double localCharge;
04704     for (atmIndx = 0; atmIndx < msg->numAllAtoms; atmIndx++) {
04705         
04706         iret = fscanf(outputFile,"%lf %lf %lf %lf\n", 
04707                       localForce+0,
04708                       localForce+1,
04709                       localForce+2,
04710                       &localCharge);
04711         if ( iret != 4 ) {
04712             NAMD_die("Error reading forces and charge from QM results file.");
04713         }
04714         
04715         // If we are reading charges and forces on QM atoms, store
04716         // them directly.
04717         if (atmIndx < msg->numQMAtoms ) {
04718             
04719             resForce[atmIndx].force.x = localForce[0];
04720             resForce[atmIndx].force.y = localForce[1];
04721             resForce[atmIndx].force.z = localForce[2];
04722             
04723             atmP[atmIndx].charge = localCharge;
04724             resForce[atmIndx].charge = localCharge;
04725         }
04726         
04727         // If we are reading charges and forces applied on Dummy atoms,
04728         // place them on the appropriate QM and MM atoms to conserve energy.
04729         
04730         // This implementation was based on the description in 
04731         // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
04732         if ( msg->numQMAtoms <= atmIndx &&
04733         atmIndx < msg->numAllAtoms ) {
04734             
04735             // We keep the calculated charge in the dummy atom
04736             // so that we can calculate electrostatic forces later on.
04737             atmP[atmIndx].charge = localCharge;
04738             
04739             // If we are reading charges from Dummy atoms,
04740             // place them on the appropriate QM atom.
04741             // The dummy atom points to the QM atom to which it is bound.
04742             int qmInd = atmP[atmIndx].bountToIndx ;
04743             resForce[qmInd].charge += localCharge;
04744             
04745             // The dummy atom points to the QM atom to which it is bound.
04746             // The QM atom and the MM atom (in a QM-MM bond) point to each other.
04747             int mmInd = atmP[qmInd].bountToIndx ;
04748             
04749             Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
04750             Real mmqmDist = dir.length() ;
04751             
04752             Real linkDist = Vector(atmP[atmIndx].position - 
04753                             atmP[qmInd].position).length() ;
04754             
04755             Force mmForce(0), qmForce(0), 
04756                 linkForce(localForce[0], localForce[1], localForce[2]);
04757             
04758             Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
04759             // Unit vectors
04760             Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
04761             Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
04762             Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
04763             Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
04764             
04765             qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv + 
04766                         (xDelta)*base) )*xuv;
04767             
04768             qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv + 
04769                         (yDelta)*base) )*yuv;
04770             
04771             qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv + 
04772                         (zDelta)*base) )*zuv;
04773             
04774             
04775             mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
04776                         (xDelta)*base) )*xuv;
04777             
04778             mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
04779                         (yDelta)*base) )*yuv;
04780             
04781             mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
04782                         (zDelta)*base) )*zuv;
04783             
04784             resForce[qmInd].force += qmForce;
04785             resForce[msg->numQMAtoms + mmInd].force += mmForce;
04786         }
04787     }
04788     
04789     // The initial point charge index for the force message is the number of QM
04790     // atoms, since the dummy atoms have no representation in NAMD
04791     int pcIndx = msg->numQMAtoms;
04792     
04793     if (usrPCnum > 0) {
04794         // We loop over point charges, reading the total electrostatic force 
04795         // applied on them by the QM region.
04796         // We redistribute the forces applied over virtual point
04797         // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
04798         for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
04799             
04800             Force totalForce(0);
04801             
04802             // No force was applied to the QM region due to this charge, since it
04803             // was skipped when writing down point charges to the QM software, so it
04804             // does not receive any back from the QM region. It must be an MM1 atom
04805             // from a QM-MM bond.
04806             if (pcP[i].type == QMPCTYPE_IGNORE)
04807                 continue;
04808             
04809             iret = fscanf(outputFile,"%lf %lf %lf\n", 
04810                            &totalForce[0], &totalForce[1], &totalForce[2]);
04811             if ( iret != 3 ) {
04812                 NAMD_die("Error reading PC forces from QM results file.");
04813             }
04814             
04815             if (pcP[i].type == QMPCTYPE_CLASSICAL) {
04816                 // If we already ignored MM1 charges, we take all other 
04817                 // non-virtual charges and apply forces directly to them.
04818                 resForce[pcIndx].force += totalForce;
04819             }
04820             else {
04821                 // If we are handling virtual PC, we distribute the force over
04822                 // MM1 and MM2.
04823                 
04824                 Force mm1Force(0), mm2Force(0);
04825                 
04826                 // Virtual PC are bound to MM2.
04827                 int mm2Indx = pcP[i].bountToIndx;
04828                 // MM2 charges are bound to MM1.
04829                 int mm1Indx = pcP[mm2Indx].bountToIndx;
04830                 
04831                 Real Cq = pcP[i].dist;
04832                 
04833                 mm1Force = (1-Cq)*totalForce ;
04834                 mm2Force = Cq*totalForce ;
04835                 
04836                 resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
04837                 resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
04838             }
04839             
04840             
04841         }
04842     }
04843     
04844     fclose(outputFile);
04845     delete [] line;
04846     
04847     // In case charges are not to be read form the QM software output,
04848     // we load the origianl atom charges.
04849     if (msg->qmAtmChrgMode == QMCHRGNONE) {
04850         int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
04851         const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
04852         Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
04853         
04854         atmIndx = 0 ;
04855         for (; atmIndx < msg->numQMAtoms; atmIndx++) {
04856             
04857             // Loops over all QM atoms (in all QM groups) comparing their global indices
04858             for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
04859                 
04860                 if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
04861                     
04862                     atmP[atmIndx].charge = qmAtmChrg[qmIter];
04863                     resForce[atmIndx].charge = qmAtmChrg[qmIter];
04864                     
04865                     break;
04866                 }
04867                 
04868             }
04869             
04870         }
04871         for (size_t j=msg->numQMAtoms; j<msg->numAllAtoms; ++j ) {
04872             atmP[j].charge = 0;
04873         }
04874     }
04875     
04876     // remove force file
04877 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
04878 //     iret = remove(outputFileName);
04879 //     if ( iret ) { NAMD_err(0); }
04880     
04881     if (usrPCnum == 0) {
04882         DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
04883         
04884         atmP = msg->data ;
04885         pcP = msg->data + msg->numAllAtoms ;
04886         
04887         // We only loop over point charges from real atoms, ignoring the ones 
04888         // created to handle QM-MM bonds.
04889         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
04890             
04891             // No force was applied to the QM region due to this charge, so it 
04892             // does not receive any back from the QM region. It must be an MM1 atom
04893             // from a QM-MM bond.
04894             if (pcP[i].type == QMPCTYPE_IGNORE)
04895                 continue;
04896             
04897             Force totalForce(0);
04898             
04899             BigReal pntCharge = pcP[i].charge;
04900             
04901             Position posMM = pcP[i].position ;
04902             
04903             for (size_t j=0; j<msg->numAllAtoms; ++j ) {
04904                 
04905                 BigReal qmCharge = atmP[j].charge ;
04906                 
04907                 BigReal force = pntCharge*qmCharge*constants ;
04908                 
04909                 Position rVec = posMM - atmP[j].position ;
04910                 
04911                 force /= rVec.length2();
04912                 
04913                 // We accumulate the total force felt by a point charge
04914                 // due to the QM charge distribution. This total force
04915                 // will be applied to the point charge if it is a real one,
04916                 // or will be distirbuted over MM1 and MM2 point charges, it 
04917                 // this is a virtual point charge.
04918                 totalForce += force*rVec.unit();
04919             }
04920             
04921             if (pcP[i].type == QMPCTYPE_CLASSICAL) {
04922                 // If we already ignored MM1 charges, we take all other 
04923                 // non-virtual charges and apply forces directly to them.
04924                 resForce[pcIndx].force += totalForce;
04925             }
04926             else {
04927                 // If we are handling virtual PC, we distribute the force over
04928                 // MM1 and MM2.
04929                 
04930                 Force mm1Force(0), mm2Force(0);
04931                 
04932                 // Virtual PC are bound to MM2.
04933                 int mm2Indx = pcP[i].bountToIndx;
04934                 // MM2 charges are bound to MM1.
04935                 int mm1Indx = pcP[mm2Indx].bountToIndx;
04936                 
04937                 Real Cq = pcP[i].dist;
04938                 
04939                 mm1Force = (1-Cq)*totalForce ;
04940                 mm2Force = Cq*totalForce ;
04941                 
04942                 resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
04943                 resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
04944             }
04945             
04946         }
04947     }
04948     
04949     // Adjusts forces from PME, canceling contributions from the QM and 
04950     // direct Coulomb forces calculated here.
04951     if (msg->PMEOn) {
04952         
04953         DebugM(1,"Correcting forces and energy for PME.\n");
04954         
04955         ewaldcof = msg->PMEEwaldCoefficient;
04956         BigReal TwoBySqrtPi = 1.12837916709551;
04957         pi_ewaldcof = TwoBySqrtPi * ewaldcof;
04958         
04959         for (size_t i=0; i < msg->numQMAtoms; i++) {
04960             
04961             BigReal p_i_charge = atmP[i].charge ;
04962             Position pos_i = atmP[i].position ;
04963             
04964             for (size_t j=i+1; j < msg->numQMAtoms; j++) {
04965                 
04966                 BigReal p_j_charge = atmP[j].charge ;
04967                 
04968                 Position pos_j = atmP[j].position ;
04969                 
04970                 BigReal r = Vector(pos_i - pos_j).length() ;
04971                 
04972                 BigReal tmp_a = r * ewaldcof;
04973                 BigReal tmp_b = erfc(tmp_a);
04974                 BigReal corr_energy = tmp_b;
04975                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
04976                 
04977 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
04978                 BigReal recip_energy = (1-tmp_b)/r;
04979                 
04980 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
04981                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
04982                 
04983                 const BigReal kq_i = p_i_charge * constants;
04984                 
04985                 // Final force and energy correction for this pair of atoms.
04986                 BigReal energy = kq_i * p_j_charge * recip_energy ;
04987                 
04988                 Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
04989                 
04990                 // The force is *subtracted* from the total force acting on
04991                 // both atoms. The sign on fixForce corrects the orientation
04992                 // of the subtracted force.
04993 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
04994 //                     << std::endl);
04995 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
04996 //                     << std::endl);
04997 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
04998                 resForce[i].force -= fixForce ;
04999                 resForce[j].force -= -1*fixForce;
05000                 
05001                 // The energy is *subtracted* from the total energy calculated here.
05002                 resMsg->energyCorr -= energy;
05003             }
05004         }
05005         
05006         pcIndx = msg->numQMAtoms;
05007         // We only loop over point charges from real atoms, ignoring the ones 
05008         // created to handle QM-MM bonds.
05009         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
05010             
05011             BigReal p_i_charge = pcPpme[i].charge;
05012             Position pos_i = pcPpme[i].position ;
05013             
05014             const BigReal kq_i = p_i_charge * constants;
05015             
05016             Force fixForce = 0;
05017             
05018             for (size_t j=0; j<msg->numQMAtoms; ++j ) {
05019                 
05020                 BigReal p_j_charge = atmP[j].charge ;
05021                 
05022                 Position pos_j = atmP[j].position ;
05023                 
05024                 BigReal r = Vector(pos_i - pos_j).length() ;
05025                 
05026                 BigReal tmp_a = r * ewaldcof;
05027                 BigReal tmp_b = erfc(tmp_a);
05028                 BigReal corr_energy = tmp_b;
05029                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
05030                 
05031 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
05032                 BigReal recip_energy = (1-tmp_b)/r;
05033                 
05034 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
05035                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
05036                 
05037                 // Final force and energy correction for this pair of atoms.
05038                 BigReal energy = kq_i * p_j_charge * recip_energy ;
05039                 
05040                 fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
05041                 
05042                 // The energy is *subtracted* from the total energy calculated here.
05043                 resMsg->energyCorr -= energy;
05044                 
05045             }
05046             
05047             // The force is *subtracted* from the total force acting on
05048                 // the point charge.
05049 //                 DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
05050 //                     << std::endl);
05051 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
05052             resForce[pcIndx].force -= kq_i*fixForce ;
05053         }
05054         
05055     }
05056     
05057     DebugM(1,"Determining virial...\n");
05058     
05059     // Calculates the virial contribution form the forces on QM atoms and 
05060     // point charges calculated here.
05061     for (size_t i=0; i < msg->numQMAtoms; i++) {
05062         // virial used by NAMD is -'ve of normal convention, so reverse it!
05063         // virial[i][j] in file should be sum of -1 * f_i * r_j
05064         for ( int k=0; k<3; ++k )
05065             for ( int l=0; l<3; ++l )
05066                 resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
05067     }
05068     
05069     pcIndx = msg->numQMAtoms; // Index in the real PC force array.
05070     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
05071         // virial used by NAMD is -'ve of normal convention, so reverse it!
05072         // virial[i][j] in file should be sum of -1 * f_i * r_j
05073         for ( int k=0; k<3; ++k )
05074             for ( int l=0; l<3; ++l )
05075                 resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
05076     }
05077     
05078     DebugM(1,"End of QM processing. Sending result message.\n");
05079     
05080     // Send message to rank zero with results.
05081     QMProxy[0].recvQMRes(resMsg);
05082     
05083     if (msg->switching && msg->PMEOn)
05084         delete [] pcPpme;
05085     
05086     delete msg;
05087     return ;
05088 }
05089 
05090 
05091 void ComputeQMMgr::pntChrgSwitching(QMGrpCalcMsg *msg, QMAtomData *pcPpme) {
05092     
05093     // We apply a switching function to the point charges so that there is a 
05094     // smooth decay of the electrostatic environment seen by the QM system.
05095     
05096     BigReal cutoff2 = msg->cutoff*msg->cutoff;
05097     BigReal swdist = msg->swdist;
05098     
05099     SortedArray<pntChrgDist> sortedDists;
05100     
05101     DebugM(1,"Initiating point charge switching and processing in rank " 
05102     << CkMyPe() << "\n" ) ;
05103     
05104     QMAtomData *pcP = msg->data + msg->numAllAtoms;
05105     
05106 #ifdef DEBUG_QM
05107     Real PCScaleCharge = 0;
05108     for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
05109         PCScaleCharge += pcP[i].charge;
05110     }
05111     DebugM(4, "The initial total Point-Charge charge is " << PCScaleCharge
05112             << " before scaling type: " << msg->switchType << "\n" );
05113 #endif
05114     
05115     // If PME is on, we return an unchanged vector of charges so that a 
05116     // PME correction can be calculated.
05117     if (msg->PMEOn) {
05118         for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
05119             pcPpme[i] = pcP[i];
05120         }
05121     }
05122     
05123     switch (msg->switchType) {
05124         
05125         case QMPCSCALESHIFT:
05126         {
05127             
05128             // We store all point charges so that only the furthest away
05129             // are changed in PC schemes "round" or "zero"
05130             for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
05131                 sortedDists.add( pntChrgDist(i, pcP[i].dist) ) ;
05132                 
05133                 // https://nmr.cit.nih.gov/xplor-nih/xplorMan/node119.html
05134                 // Applying X-PLOR shifting formula:
05135                 // F = Qi*Qj*(C/(epsilon*r))*(1 - r^2/cutoff^2)^2
05136                 Real dist2 = pcP[i].dist*pcP[i].dist ;
05137                 dist2 /= cutoff2 ;
05138                 Real coef = 1- dist2;
05139                 pcP[i].charge *= coef*coef ;
05140             }
05141             
05142         } break;
05143         
05144         case QMPCSCALESWITCH:
05145         {
05146             for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
05147                 
05148                 // We store the point charges which are beiond the switching
05149                 // distance.
05150                 if (pcP[i].dist > swdist) {
05151                     sortedDists.add( pntChrgDist(i, pcP[i].dist) ) ;
05152                     
05153                     // https://nmr.cit.nih.gov/xplor-nih/xplorMan/node118.html
05154                     // Applying the XPLOR switching formula:
05155                     // (r^2 - cutoff^2)^2*(cutoff^2 + 2*r^2 - 3*swdits^2)/(cutoff^2 - swdits^2)^3
05156                     Real dist2 = pcP[i].dist*pcP[i].dist ;
05157                     Real swdist2 = swdist*swdist;
05158                     Real coef = (dist2 - cutoff2)*(dist2 - cutoff2) ;
05159                     coef *= (cutoff2 + 2*dist2 - 3*swdist2) ;
05160                     coef /= (cutoff2 - swdist2)*(cutoff2 - swdist2)*(cutoff2 - swdist2);
05161                     pcP[i].charge *= coef ;
05162                 }
05163             }
05164         } break;
05165     }
05166     
05167 #ifdef DEBUG_QM
05168     PCScaleCharge = 0;
05169     for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
05170         PCScaleCharge += pcP[i].charge;
05171     }
05172     DebugM(4, "The final total Point-Charge charge is " << PCScaleCharge
05173             << " after scaling.\n" );
05174 #endif
05175     
05176     DebugM(4, sortedDists.size() 
05177     << " point charges were selected for point charge scheme " << msg->pcScheme << "\n" );
05178     
05179     Real totalPCCharge = 0, correction = 0;
05180     switch (msg->pcScheme) {
05181         
05182         case QMPCSCHEMEROUND:
05183         {
05184             for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
05185                 totalPCCharge += pcP[i].charge;
05186             }
05187             DebugM(4, "The total Point-Charge charge is " << totalPCCharge 
05188             << "\n" );
05189             
05190             if ((fabsf(roundf(totalPCCharge) - totalPCCharge) <= 0.001f) ) {
05191                 DebugM(4, "Charge is already a whole number!\n" );
05192             } 
05193             else {
05194                 correction = roundf(totalPCCharge) -totalPCCharge ;
05195                 DebugM(4, "Adding to system the charge: " << correction << "\n" );
05196             }
05197         } break;
05198         
05199         case QMPCSCHEMEZERO:
05200         {
05201             for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
05202                 totalPCCharge += pcP[i].charge;
05203             }
05204             DebugM(4, "The total Point-Charge charge is " << totalPCCharge << "\n");
05205             
05206             DebugM(4, "Total QM system charge is: " << msg->charge << "\n" );
05207             
05208             correction = -1*(totalPCCharge + msg->charge);
05209             if ((fabsf(correction) <= 0.001f) ) {
05210                 correction = 0;
05211                 DebugM(4, "Total QM + PC charge is already zero!\n" );
05212             }
05213             else
05214                 DebugM(4, "Adding a charge of " << correction << " to the system\n");
05215             
05216         } break;
05217     }
05218     
05219     if (correction != 0) {
05220         
05221         int maxi = sortedDists.size(), mini = sortedDists.size()/2;
05222         Real fraction = correction/(maxi - mini); 
05223         
05224         DebugM(4, "The charge fraction added to the " << (maxi - mini)
05225         << " most distant point charges is " << fraction << "\n");
05226         
05227         for (size_t i=mini; i<maxi ; i++) {
05228             
05229             pcP[sortedDists[i].index].charge += fraction ;
05230             
05231         }
05232         
05233         #ifdef DEBUG_QM
05234         totalPCCharge = 0;
05235         for ( size_t i=0; i<msg->numRealPntChrgs; i++ ) {
05236             totalPCCharge += pcP[i].charge;
05237         }
05238         DebugM(4, "The total Point-Charge charge is " << totalPCCharge 
05239         << "\n");
05240         #endif
05241     }
05242     
05243 }
05244 
05245 void ComputeQMMgr::lssPrepare() {
05246     
05247     lssTotRes = 0;
05248     lssResMass = 0;
05249     
05250     DebugM (4, "Preparing LSS for " << numQMGrps << " QM groups.\n" )
05251     
05252     for (int i=0; i<numQMGrps; i++) {
05253         lssTotRes += qmLSSSize[i];
05254     }
05255     
05256     lssPos = new Position[lssTotRes];
05257     
05258     grpIDResNum.resize(numQMGrps);
05259     
05260     if (simParams->qmLSSMode == QMLSSMODECOM) {
05261         
05262         lssGrpRefMass.resize(numQMGrps);
05263         
05264         for (int i=0; i<qmLSSResSize; i++)
05265             lssResMass += qmLSSMass[i];
05266         
05267         DebugM(4, "Total residue mass is " << lssResMass << "\n" )
05268     }
05269     
05270     // Get all atom IDs of solvent QM atoms, per group.
05271     int solvResBeg = 0, refAtmBeg = 0, locIter = 0, solvIndx = 0;
05272     int *lssIndxs, *refAtmsIndxs ;
05273     Mass *lssMasses, *refAtmMasses;
05274     while (locIter < numQMGrps) {
05275         lssIndxs = qmLSSIdxs + solvResBeg;
05276         
05277         DebugM (4, "Loading atom IDs for QM group " << locIter 
05278             << " with " << qmLSSSize[locIter]
05279             << " solvent molecules.\n" )
05280         
05281         
05282         switch (simParams->qmLSSMode) {
05283         
05284         case QMLSSMODECOM:
05285         {
05286             lssMasses = qmLSSMass + solvResBeg;
05287             
05288             // Loads data on QM solvent residues and their atoms.
05289             for (int i=0; i<qmLSSSize[locIter]; i++) {
05290                 
05291                 lssPos[solvIndx] = 0;
05292                 
05293                 Debug( iout << "Getting atom IDs from QM solvent molecule " 
05294                 << solvIndx << "\n") ;
05295                 for (int j=0; j<qmLSSResSize; j++) {
05296                     
05297                     int atmID = lssIndxs[i*qmLSSResSize + j];
05298                     Mass atmMass = lssMasses[i*qmLSSResSize + j];
05299                     Debug( iout << atmID << " (" << atmMass << ") ") ;
05300                     
05301                     grpIDResNum[locIter].insert(atmLSSData(atmID, LSSDataStr(solvIndx,atmMass)));
05302                     
05303                 }
05304                 solvIndx++;
05305                 Debug( iout << "\n" << endi );
05306             }
05307             
05308             // Loads data on the mass of QM atoms which will be used to calculate 
05309             // the COM for solvent selection.
05310             refAtmsIndxs = qmLSSRefIDs + refAtmBeg;
05311             refAtmMasses = molPtr->get_qmLSSRefMass() + refAtmBeg;
05312             for (int i=0; i<qmLSSRefSize[locIter]; i++) {
05313                 lssGrpRefMass[locIter].insert( 
05314                     refLSSData( refAtmsIndxs[i], refAtmMasses[i] ) 
05315                 ) ;
05316             }
05317             refAtmBeg += qmLSSRefSize[locIter] ;
05318         } break ;
05319         
05320         case QMLSSMODEDIST:
05321         {
05322             // Loads data on QM solvent residues and their atoms.
05323             for (int i=0; i<qmLSSSize[locIter]; i++) {
05324                 
05325                 lssPos[solvIndx] = 0;
05326                 
05327                 Debug( iout << "Getting atom IDs from QM solvent molecule " 
05328                 << solvIndx << "\n") ;
05329                 for (int j=0; j<qmLSSResSize; j++) {
05330                     
05331                     int atmID = lssIndxs[i*qmLSSResSize + j];
05332                     Debug( iout << atmID << " ") ;
05333                     
05334                     grpIDResNum[locIter].insert(atmLSSData(atmID, LSSDataStr(solvIndx,0)));
05335                     
05336                 }
05337                 solvIndx++;
05338                 Debug( iout << "\n" << endi );
05339             }
05340             
05341         } break ;
05342         
05343         }
05344         
05345         solvResBeg += qmLSSSize[locIter]*qmLSSResSize ;
05346         locIter++;
05347     }
05348     
05349     return ;
05350 }
05351 
05352 void ComputeQMMgr::lssUpdate(int grpIter, QMAtmVec& grpQMAtmVec, 
05353                              QMPCVec& grpPntChrgVec) {
05354     
05355     SortedArray<lssDistSort> solvDist;
05356     
05357     Position refCOM(0) ;
05358     Mass totMass = 0;
05359     
05360     DebugM(3, "LSS UPDATE...\n")
05361     
05362     int solvResBeg = 0 ;
05363     for (int i=0; i<grpIter; i++)
05364         solvResBeg += qmLSSSize[i] ;
05365     
05366     switch (simParams->qmLSSMode ) {
05367         
05368     case QMLSSMODECOM:
05369     {
05370         DebugM(3, "Using COM for LSS in group " << grpIter << "\n")
05371         
05372         // Determined the reference center of mass for this group.
05373         for(int i=0; i<grpQMAtmVec.size(); i++) {
05374             
05375             auto it = lssGrpRefMass[grpIter].find(grpQMAtmVec[i].id);
05376             if ( it != lssGrpRefMass[grpIter].end() ) {
05377                 refCOM += grpQMAtmVec[i].position*it->second ;
05378                 totMass += it->second ;
05379             }
05380         }
05381         
05382         refCOM /= totMass;
05383         DebugM ( 3, "Reference COM position: " << refCOM << "\n");
05384         
05385         // Initialize the positions of all COM of quantum residues 
05386         for (int solvIter=solvResBeg; solvIter<solvResBeg+qmLSSSize[grpIter]; solvIter++) {
05387             lssPos[solvIter] = 0 ;
05388         }
05389         
05390         DebugM(3, "Calculating distance of QM solvent COM from reference COM of group.\n")
05391         
05392         // Temporary handler of lssDistSort structures, while we accumulate atoms
05393         // on their respective residues and calculate distances.
05394         std::map<int,lssDistSort > resQMDist ;
05395         
05396         // Initiates COM determination for all QM solvent molecules.
05397         for(int i=0; i<grpQMAtmVec.size(); i++) {
05398             auto it = grpIDResNum[grpIter].find(grpQMAtmVec[i].id) ;
05399             if (it != grpIDResNum[grpIter].end()) {
05400                 lssPos[it->second.resIndx] += grpQMAtmVec[i].position*it->second.mass ;
05401                 
05402                 // tries to find a residue number
05403                 auto itRes = resQMDist.find(it->second.resIndx) ;
05404                 if (itRes == resQMDist.end() ) {
05405                     resQMDist.insert(std::pair<int,lssDistSort >(
05406                         it->second.resIndx, lssDistSort(QMLSSQMRES, -1) ));
05407                 }
05408                 
05409                 // For each classical residue ID, we compile a list of atom IDs which 
05410                 // are found among point charges
05411 //                 resQMDist[it->second.resIndx].idVec.push_back(grpQMAtmVec[i].id) ;
05412 //                 resQMDist[it->second.resIndx].indxVec.push_back(i) ;
05413                 
05414                 resQMDist[it->second.resIndx].idIndx.add(idIndxStr(grpQMAtmVec[i].id,i)) ;
05415             }
05416         }
05417         
05418         DebugM(3, "QM Solvent molecules " << solvResBeg 
05419         << " to " << solvResBeg+qmLSSSize[grpIter] << "\n")
05420         
05421         
05422         for (int solvIter=solvResBeg; solvIter<solvResBeg+qmLSSSize[grpIter]; solvIter++) {
05423             Real dist = Vector((lssPos[solvIter]/lssResMass) - refCOM).length() ;
05424             resQMDist[solvIter].dist = dist ;
05425             solvDist.add(resQMDist[solvIter]);
05426         }
05427         
05428         #ifdef DEBUG_QM
05429         DebugM(3, "We loaded the following QM solvent residues and distances:\n")
05430         for (int i=0; i<solvDist.size(); i++) {
05431             iout << i << ") type: " << solvDist[i].type 
05432             << " dist " << solvDist[i].dist
05433             << " IDs: " ;
05434             for (int j=0; j<solvDist[i].idIndx.size(); j++) 
05435                 iout << solvDist[i].idIndx[j].ID << " " ;
05436             iout << "\n" << endi;
05437         }
05438         #endif
05439         
05440         // Get Center Of Mass distance for all (whole) solvent molecules among 
05441         // *point charges*, comparing against the non-solvent QM atoms.
05442         
05443         // This will count how many PCs are associated with each solvent residue,
05444         // and associate their atom ID with that residue.
05445         std::map<int,lssDistSort > resPCSize ;
05446         
05447         for(int i=0; i<grpPntChrgVec.size(); i++) {
05448             
05449             // This maps aomt IDs with a global classical residue ID.
05450             auto it = molPtr->get_qmMMSolv().find(grpPntChrgVec[i].id) ;
05451             if (it != molPtr->get_qmMMSolv().end()) {
05452                 
05453                 // tries to find a residue number
05454                 auto itRes = resPCSize.find(it->second) ;
05455                 if (itRes == resPCSize.end() ) {
05456                     resPCSize.insert(std::pair<int,lssDistSort >(
05457                         it->second, lssDistSort(QMLSSCLASSICALRES, -1) ));
05458                 }
05459                 
05460                 // For each classical residue ID, we compile a list of atom IDs which 
05461                 // are found among point charges
05462 //                 resPCSize[it->second].idVec.push_back(grpPntChrgVec[i].id) ;
05463 //                 resPCSize[it->second].indxVec.push_back(i) ;
05464                 
05465                 resPCSize[it->second].idIndx.add(idIndxStr(grpPntChrgVec[i].id,i)) ;
05466             }
05467         }
05468         
05469         // Now we check which classical solvent molecules are complete,
05470         // and compute the COM for the complete ones, while ignoring the 
05471         // incomplete ones.
05472         for (auto it=resPCSize.begin(); it!=resPCSize.end(); it++) {
05473             
05474             if (it->second.idIndx.size() == qmLSSResSize) {
05475                 
05476                 Position currCOM(0);
05477                 Mass totalMass = 0;
05478                 
05479     //             iout << "Found complete classical residue " << it->first << "\n" << endi;
05480                 
05481                 for (int i=0; i<it->second.idIndx.size(); i++) {
05482     //                 iout << "AtomID " << it->second.idIndx[i] .ID
05483     //                 << " indxVec " << it->second.idIndx[i].indx
05484     //                 << " mass " << grpPntChrgVec[it->second.idIndx[i].indx].mass
05485     //                 << " position " << grpPntChrgVec[it->second.idIndx[i].indx].position << "\n" << endi;
05486                     currCOM += grpPntChrgVec[it->second.idIndx[i].indx].position*grpPntChrgVec[it->second.idIndx[i].indx].mass;
05487                     totalMass += grpPntChrgVec[it->second.idIndx[i].indx].mass;
05488                 }
05489                 currCOM /= totalMass;
05490                 
05491                 Real currDist = Vector(currCOM - refCOM).length() ;
05492                 
05493                 it->second.dist = currDist ;
05494                 
05495                 solvDist.add(it->second) ;
05496             }
05497             
05498         }
05499     } break ;
05500     
05501     case QMLSSMODEDIST:
05502     {
05503         DebugM(3, "Using minimal distances for LSS in group " << grpIter << "\n")
05504         
05505         DebugM(3, "QM Solvent molecules " << solvResBeg 
05506         << " to " << solvResBeg+qmLSSSize[grpIter] << "\n")
05507         
05508         // List of QM indices which will be used as reference for distance calculation.
05509         ResizeArray<int> qmRefIndx ;
05510         
05511         // Temporary handler of lssDistSort structures, while we accumulate atoms
05512         // on their respective residues and calculate distances.
05513         std::map<int,lssDistSort > resQMDist ;
05514         
05515         // Initiates COM determination for all QM solvent molecules.
05516         for(int i=0; i<grpQMAtmVec.size(); i++) {
05517             
05518             auto it = grpIDResNum[grpIter].find(grpQMAtmVec[i].id) ;
05519             if (it != grpIDResNum[grpIter].end()) {
05520                 
05521                 // tries to find a residue number
05522                 auto itRes = resQMDist.find(it->second.resIndx) ;
05523                 if (itRes == resQMDist.end() ) {
05524                     resQMDist.insert(std::pair<int,lssDistSort >(
05525                         it->second.resIndx, lssDistSort(QMLSSQMRES, -1) ));
05526                 }
05527                 
05528                 // For each classical residue ID, we compile a list of atom IDs which 
05529                 // are found among point charges
05530 //                 resQMDist[it->second.resIndx].idVec.push_back(grpQMAtmVec[i].id) ;
05531 //                 resQMDist[it->second.resIndx].indxVec.push_back(i) ;
05532                 
05533                 resQMDist[it->second.resIndx].idIndx.add(idIndxStr(grpQMAtmVec[i].id,i)) ;
05534                 
05535             }
05536             else {
05537                 qmRefIndx.add(i) ;
05538             }
05539         }
05540         
05541         // We now calculate the shortest distance between a QM solvent residue 
05542         // and any non-solvent QM atom.
05543         for (auto it=resQMDist.begin(); it != resQMDist.end(); it++) {
05544             
05545             // We prime the residue distance with the first non-solvent
05546             // QM atom.
05547             it->second.dist = Vector(
05548                     grpQMAtmVec[it->second.idIndx[0].indx].position - 
05549                     grpQMAtmVec[qmRefIndx[0]].position
05550                     ).length() ;
05551             
05552             for (int i=0; i<it->second.idIndx.size(); i++) {
05553                 
05554                 for(int j=0; j<qmRefIndx.size(); j++) {
05555                     Real currDist = Vector(
05556                         grpQMAtmVec[it->second.idIndx[i].indx].position - 
05557                         grpQMAtmVec[qmRefIndx[j]].position
05558                         ).length() ;
05559                     
05560                     if (currDist < it->second.dist)
05561                         it->second.dist = currDist;
05562                 }
05563             }
05564             
05565             // After determining the distance of this QM solvent residue,
05566             // we add it to the sorted list.
05567             solvDist.add(it->second) ;
05568         }
05569         
05570                 
05571         #ifdef DEBUG_QM
05572         DebugM(3, "We loaded the following QM solvent residues and distances:\n")
05573         for (int i=0; i<solvDist.size(); i++) {
05574             iout << i << ") type: " << solvDist[i].type 
05575             << " dist " << solvDist[i].dist
05576             << " IDs: " ;
05577             for (int j=0; j<solvDist[i].idIndx.size(); j++) 
05578                 iout << solvDist[i].idIndx[j].ID << " " ;
05579             iout << "\n" << endi;
05580         }
05581         #endif
05582         
05583         // Get shortest distance for all (whole) solvent molecules among 
05584         // *point charges*, comparing against the non-solvent QM atoms.
05585         
05586         // This will count how many PCs are associated with each solvent residue,
05587         // and associate their atom ID with that residue.
05588         std::map<int,lssDistSort > resPCSize ;
05589         
05590         for(int i=0; i<grpPntChrgVec.size(); i++) {
05591             
05592             // This maps aomt IDs with a global classical residue ID.
05593             auto it = molPtr->get_qmMMSolv().find(grpPntChrgVec[i].id) ;
05594             if (it != molPtr->get_qmMMSolv().end()) {
05595                 
05596                 // tries to find a residue number
05597                 auto itRes = resPCSize.find(it->second) ;
05598                 if (itRes == resPCSize.end() ) {
05599                     resPCSize.insert(std::pair<int,lssDistSort >(
05600                         it->second, lssDistSort(QMLSSCLASSICALRES, -1) ));
05601                 }
05602                 
05603                 // For each classical residue ID, we compile a list of atom IDs which 
05604                 // are found among point charges
05605 //                 resPCSize[it->second].idVec.push_back(grpPntChrgVec[i].id) ;
05606 //                 resPCSize[it->second].indxVec.push_back(i) ;
05607                 
05608                 resPCSize[it->second].idIndx.add(idIndxStr(grpPntChrgVec[i].id,i)) ;
05609             }
05610         }
05611         
05612         // Now we check which classical solvent molecules are complete,
05613         // and compute the COM for the complete ones, while ignoring the 
05614         // incomplete ones.
05615         for (auto it=resPCSize.begin(); it!=resPCSize.end(); it++) {
05616             
05617             if (it->second.idIndx.size() == qmLSSResSize) {
05618                 
05619                 // We prime the residue distance with the first non-solvent
05620                 // QM atom.
05621                 it->second.dist = Vector(
05622                         grpPntChrgVec[it->second.idIndx[0].indx].position - 
05623                         grpQMAtmVec[qmRefIndx[0]].position
05624                         ).length() ;
05625                 
05626                 for (int i=0; i<it->second.idIndx.size(); i++) {
05627                     
05628                     for(int j=0; j<qmRefIndx.size(); j++) {
05629                         Real currDist = Vector(
05630                             grpPntChrgVec[it->second.idIndx[i].indx].position - 
05631                             grpQMAtmVec[qmRefIndx[j]].position
05632                             ).length() ;
05633                         
05634                         if (currDist < it->second.dist)
05635                             it->second.dist = currDist;
05636                     }
05637                 }
05638                 
05639                 solvDist.add(it->second) ;
05640             }
05641             
05642         }
05643         
05644     } break ;
05645     
05646     } // End switch
05647     
05648     #ifdef DEBUG_QM
05649     DebugM(3, "Final selection of solvent residues and distances:\n")
05650     for (int i=0; i<qmLSSSize[grpIter]; i++) {
05651         std::string typeS ;
05652         if (solvDist[i].type != QMLSSCLASSICALRES ) 
05653             typeS = "QM" ;
05654         else 
05655             typeS = "Classical";
05656         iout << i << ") type: " << typeS
05657         << " dist " << solvDist[i].dist
05658         << " IDs: " ;
05659         for (int j=0; j<solvDist[i].idIndx.size(); j++) 
05660             iout << solvDist[i].idIndx[j].ID << " " ;
05661         iout << "\n" << endi;
05662     }
05663     #endif
05664     
05665     // Compare COM distances of QM and Classical solvent molecules, creating a
05666     // substitution list, atmID by atmID.
05667     
05668     DebugM(3, "Determining residues to be swaped...\n")
05669     
05670     ResizeArray<lssDistSort> nearMM, farQM ;
05671     
05672     for (int resIter = 0; resIter < qmLSSSize[grpIter] ; resIter++) {
05673         if (solvDist[resIter].type == QMLSSCLASSICALRES) {
05674             nearMM.add(solvDist[resIter]);
05675         }
05676     }
05677     
05678     for (int resIter=qmLSSSize[grpIter]; resIter<solvDist.size(); resIter++) {
05679         if (solvDist[resIter].type == QMLSSQMRES) {
05680             farQM.add(solvDist[resIter]);
05681         }
05682         
05683         if (farQM.size() == nearMM.size()) break;
05684     }
05685     
05686     if (farQM.size() != nearMM.size())
05687         NAMD_die("Could not find complementing residues to be swapped in LSS.") ;
05688     
05689     #ifdef DEBUG_QM
05690     DebugM(3, "Removing the following QM residues:\n")
05691     for (int i=0; i<farQM.size();i++) {
05692         std::string typeS ;
05693         if (farQM[i].type != QMLSSCLASSICALRES ) 
05694             typeS = "QM" ;
05695         else 
05696             typeS = "Classical";
05697         iout << i << ") type: " << typeS
05698         << " dist " << farQM[i].dist
05699         << " IDs: " ;
05700         for (int j=0; j<farQM[i].idIndx.size(); j++) 
05701             iout << farQM[i].idIndx[j].ID << " " ;
05702         iout << "\n" << endi;
05703     }
05704     
05705     DebugM(3, "Replacing with the following Classical residues:\n")
05706     for (int i=0; i<nearMM.size();i++) {
05707         std::string typeS ;
05708         if (nearMM[i].type != QMLSSCLASSICALRES ) 
05709             typeS = "QM" ;
05710         else 
05711             typeS = "Classical";
05712         iout << i << ") type: " << typeS
05713         << " dist " << nearMM[i].dist
05714         << " IDs: " ;
05715         for (int j=0; j<nearMM[i].idIndx.size(); j++) 
05716             iout << nearMM[i].idIndx[j].ID << " " ;
05717         iout << "\n" << endi;
05718     }
05719     
05720     DebugM(3, "Building substitution array...\n")
05721     #endif
05722     
05723     iout << iINFO << "LSS is swapping " << farQM.size() << " solvent residues in QM group "
05724     << grpIter << "\n" << endi;
05725     
05726     // Now we build the array which will be sent to all nodes with force results from
05727     // this step.
05728     // Atom reassignment will be done in the next step, and will use this data.
05729     for (int i=0; i<farQM.size();i++) {
05730         
05731         for(int j=0; j<qmLSSResSize; j++) {
05732             
05733             int qIndx= farQM[i].idIndx[j].indx;
05734             int mIndx= nearMM[i].idIndx[j].indx;
05735             
05736             subsArray.add( LSSSubsDat( grpQMAtmVec[qIndx].id,
05737                                        grpPntChrgVec[mIndx].id,
05738                                        grpPntChrgVec[mIndx].vdwType,
05739                                        grpPntChrgVec[mIndx].charge ) );
05740             
05741             subsArray.add( LSSSubsDat( grpPntChrgVec[mIndx].id,
05742                                        grpQMAtmVec[qIndx].id,
05743                                        grpQMAtmVec[qIndx].vdwType,
05744                                        grpQMAtmVec[qIndx].charge ) );
05745         }
05746         
05747     }
05748     
05749     #ifdef DEBUG_QM
05750     for(int i=0; i<subsArray.size() ;i++) {
05751         DebugM(3, CkMyPe() << ") Storing LSS atom " << subsArray[i].origID
05752         << " Which will become " << subsArray[i].newID
05753         << " - " << subsArray[i].newVdWType
05754         << " - " << subsArray[i].newCharge << "\n" ); 
05755     }
05756     #endif
05757     
05758     return;
05759 }
05760 
05761 void ComputeQMMgr::calcCSMD(int grpIndx, int numQMAtoms, const QMAtomData *atmP, Force *cSMDForces) {
05762     
05763     // For each Conditional SMD instance, we update the guide particle
05764     // and calculate and apply the respective force.
05765     for ( int i=0; i < cSMDindxLen[grpIndx]; i++ ) {
05766         
05767         // Defite target distance
05768         Real tdist;
05769         Force cSMDforce(0);
05770         
05771         // cSMD instance index
05772         int cSMDit = cSMDindex[grpIndx][i] ;
05773         AtomID src = -1, trgt = -1;
05774         
05775         // Finds local indices for the atoms involved in this cSMD instance
05776         for (size_t indx=0; indx < numQMAtoms; ++indx) {
05777             
05778             if ( cSMDpairs[cSMDit][0] == atmP[indx].id )
05779                 src = indx ;
05780             
05781             if ( cSMDpairs[cSMDit][1] == atmP[indx].id )
05782                 trgt = indx ;
05783             
05784             if (src != -1 && trgt != -1)
05785                 break;
05786         }
05787         
05788         // Finds the unit vector in the direction of the target
05789         Vector trgtDir = atmP[trgt].position - atmP[src].position ;
05790         tdist = trgtDir.length();
05791         trgtDir /= tdist;
05792         
05793         // If this is the first time step, set up the guide particle.
05794         if (cSMDInitGuides < cSMDnumInstances) {
05795             cSMDguides[cSMDit] = atmP[src].position;
05796             cSMDInitGuides++;
05797         }
05798         
05799         // If the target is further away than "cutoff", advance the
05800         // guide particle and apply force.
05801         if (tdist > cSMDcoffs[cSMDit]) {
05802             
05803             cSMDguides[cSMDit] += trgtDir*cSMDVels[cSMDit];
05804             
05805             Vector guideDir = cSMDguides[cSMDit] - atmP[src].position ;
05806             
05807             // Calculate force in the direction of the guide particle.
05808             guideDir *= cSMDKs[cSMDit] ;
05809             cSMDforce = guideDir ;
05810              
05811         } 
05812         else {
05813             // If we pulled the Src towards Trgt closer than Cutoff,
05814             // we reset the guide particle to the position of Src.
05815             cSMDguides[cSMDit] = atmP[src].position;
05816         }
05817         
05818         DebugM(4, cSMDit << ") Center at  " << cSMDguides[cSMDit] << " for target distance " << 
05819         tdist << " and direction " << trgtDir << 
05820         "." << std::endl);
05821         
05822         // We now save the force in a vector of length "cSMDnumInstances".
05823         cSMDForces[cSMDit] = cSMDforce;
05824         
05825         iout << iINFO << "Calculated cSMD force vector " << cSMDforce
05826         << " (atom "  << cSMDpairs[cSMDit][0] << " to " << cSMDpairs[cSMDit][1] << ")\n" << endi;
05827     }
05828 }
05829 
05830 
05831 
05832 #ifdef DEBUG_QM
05833 
05834 void ComputeQMMgr::Write_PDB(std::string Filename, const QMGrpCalcMsg *dataMsg)
05835 {
05836     std::ofstream OutputTmpPDB ;
05837     
05838     try
05839     {
05840         
05841         OutputTmpPDB.open( Filename.c_str(), std::fstream::out );
05842         
05843         OutputTmpPDB << "REMARK Information used by NAMD to create the QM simulation." << std::endl ;
05844         OutputTmpPDB << "REMARK Occupancy: bount to index; Beta: System ID; " << std::endl ;
05845         
05846         const QMAtomData *dataP = dataMsg->data;
05847         
05848         for ( int i=0 ; i < dataMsg->numAllAtoms + dataMsg->numAllPntChrgs ; i++ )
05849         {
05850             // PDB format: http://www.wwpdb.org/documentation/format33/sect9.html
05851             
05852             // Exception: The field resName was changed from the coluns 18-20 to the columns 18-21
05853             // This allows the usage of protonated amino acids and other molecules in the PDB file.
05854             
05855             // Exception2: The field extraInfo was changed from the coluns 67 - 76 to the columns 67 - 72
05856             // This allows the usage of segments in PDB files, necessary for the propper usage of NAMD.
05857             
05858             std::string name(" x  ");
05859             std::string resName ("  uk");
05860             std::string chainName("X");
05861             std::string element("") ;
05862             if (i < dataMsg->numAllAtoms ) {
05863                 name  = dataP[i].element;
05864                 chainName = "q" ;
05865                 element = dataP[i].element;
05866                 if (dataP[i].type == QMATOMTYPE_QM)
05867                     resName = " qm " ;
05868                 else if (dataP[i].type == QMATOMTYPE_DUMMY)
05869                     resName = " dm " ;
05870             }
05871             else {
05872                 chainName = "c" ;
05873                 if (dataP[i].type == QMPCTYPE_CLASSICAL)
05874                     resName = " pc ";
05875                 else if (dataP[i].type == QMPCTYPE_IGNORE)
05876                     resName = "ipc ";
05877                 else if (dataP[i].type == QMPCTYPE_EXTRA)
05878                     resName = "vpc ";
05879             }
05880             
05881             OutputTmpPDB << "ATOM  " ; // ATOM  1 -  6
05882             
05883             OutputTmpPDB.width(5) ; // serial  7 - 11
05884             OutputTmpPDB.right ;
05885             OutputTmpPDB << i ;
05886             
05887             OutputTmpPDB << " " ; // Spacing
05888             
05889             OutputTmpPDB.width(4) ; // name  13 - 16
05890             OutputTmpPDB << name ;
05891             
05892             OutputTmpPDB << " " ; // altLoc  17
05893             
05894             OutputTmpPDB.width(4) ; // resName  18 - 21
05895             OutputTmpPDB << resName ;
05896             
05897             OutputTmpPDB.width(1) ; // chainID  22
05898             OutputTmpPDB << chainName ;
05899             
05900             OutputTmpPDB.width(4) ; // Residue Index  23 - 26
05901             OutputTmpPDB << i ;
05902             
05903             OutputTmpPDB << " " ; // iCode  27
05904             
05905             OutputTmpPDB << "   " ; // Spacing
05906             
05907             OutputTmpPDB.width(8) ; // x  31 - 38
05908             OutputTmpPDB.right ;
05909             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
05910             OutputTmpPDB.precision(3) ;
05911             OutputTmpPDB << dataP[i].position.x ;
05912             
05913             OutputTmpPDB.width(8) ; // y  39 - 46
05914             OutputTmpPDB.right ;
05915             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
05916             OutputTmpPDB.precision(3) ;
05917             OutputTmpPDB << dataP[i].position.y ;
05918             
05919             OutputTmpPDB.width(8) ; // z  47 - 54
05920             OutputTmpPDB.right ;
05921             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
05922             OutputTmpPDB.precision(3) ;
05923             OutputTmpPDB << dataP[i].position.z ;
05924             
05925             OutputTmpPDB.width(6) ; // occupancy 55 - 60
05926             OutputTmpPDB.right ;
05927             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
05928             OutputTmpPDB.precision(2) ;
05929             OutputTmpPDB << dataP[i].bountToIndx ;
05930             
05931             OutputTmpPDB.width(6) ; // tempFactor/Beta 61 - 66
05932             OutputTmpPDB.right ;
05933             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
05934             OutputTmpPDB.precision(2) ;
05935             OutputTmpPDB << dataP[i].id ;
05936             
05937             OutputTmpPDB.width(6) ; // extra information not originaly on PDB format  67 - 72
05938             OutputTmpPDB << "      " ;
05939             
05940             OutputTmpPDB.width(4) ; // segment information from NAMD, not originaly on PDB format  72 - 76
05941             OutputTmpPDB.left ;
05942             OutputTmpPDB << "QM  ";
05943             
05944             OutputTmpPDB.width(2) ; // element 77 - 78
05945             OutputTmpPDB.right ;
05946             OutputTmpPDB << element ;
05947             
05948             OutputTmpPDB.width(2) ; // charge 77 - 78
05949             OutputTmpPDB.right ;
05950             OutputTmpPDB << dataP[i].charge ;
05951             
05952             OutputTmpPDB << std::endl;
05953             
05954         }
05955         
05956         OutputTmpPDB << "END" << std::endl;
05957         
05958         OutputTmpPDB.close();
05959     }
05960     catch (...)
05961     {
05962         iout << iERROR << "Generic exception at QM write PBD." << endi ;
05963         NAMD_die("PDB write error");
05964         throw "Generic exception!" ;
05965     }
05966     return ;
05967 }
05968 
05969 void ComputeQMMgr::Write_PDB(std::string Filename, const QMCoordMsg *dataMsg)
05970 {
05971     std::ofstream OutputTmpPDB ;
05972     
05973     try
05974     {
05975         
05976         OutputTmpPDB.open( Filename.c_str(), std::fstream::out );
05977         
05978         OutputTmpPDB << "REMARK Information used by NAMD to create the QM simulation." << std::endl ;
05979         OutputTmpPDB << "REMARK Occupancy: bount to index; Beta: System ID; " << std::endl ;
05980         
05981         const ComputeQMAtom *dataP = dataMsg->coord;
05982         
05983         for ( int i=0 ; i < dataMsg->numAtoms; i++ )
05984         {
05985             // PDB format: http://www.wwpdb.org/documentation/format33/sect9.html
05986             
05987             // Exception: The field resName was changed from the coluns 18-20 to the columns 18-21
05988             // This allows the usage of protonated amino acids and other molecules in the PDB file.
05989             
05990             // Exception2: The field extraInfo was changed from the coluns 67 - 76 to the columns 67 - 72
05991             // This allows the usage of segments in PDB files, necessary for the propper usage of NAMD.
05992             
05993             std::string name(" x  ");
05994             std::string resName (" atm");
05995             std::string chainName("X");
05996             std::string element("") ;
05997             
05998             OutputTmpPDB << "ATOM  " ; // ATOM  1 -  6
05999             
06000             OutputTmpPDB.width(5) ; // serial  7 - 11
06001             OutputTmpPDB.right ;
06002             OutputTmpPDB << i ;
06003             
06004             OutputTmpPDB << " " ; // Spacing
06005             
06006             OutputTmpPDB.width(4) ; // name  13 - 16
06007             OutputTmpPDB << name ;
06008             
06009             OutputTmpPDB << " " ; // altLoc  17
06010             
06011             OutputTmpPDB.width(4) ; // resName  18 - 21
06012             OutputTmpPDB << resName ;
06013             
06014             OutputTmpPDB.width(1) ; // chainID  22
06015             OutputTmpPDB << chainName ;
06016             
06017             OutputTmpPDB.width(4) ; // Residue Index  23 - 26
06018             OutputTmpPDB << i ;
06019             
06020             OutputTmpPDB << " " ; // iCode  27
06021             
06022             OutputTmpPDB << "   " ; // Spacing
06023             
06024             OutputTmpPDB.width(8) ; // x  31 - 38
06025             OutputTmpPDB.right ;
06026             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
06027             OutputTmpPDB.precision(3) ;
06028             OutputTmpPDB << dataP[i].position.x ;
06029             
06030             OutputTmpPDB.width(8) ; // y  39 - 46
06031             OutputTmpPDB.right ;
06032             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
06033             OutputTmpPDB.precision(3) ;
06034             OutputTmpPDB << dataP[i].position.y ;
06035             
06036             OutputTmpPDB.width(8) ; // z  47 - 54
06037             OutputTmpPDB.right ;
06038             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
06039             OutputTmpPDB.precision(3) ;
06040             OutputTmpPDB << dataP[i].position.z ;
06041             
06042             OutputTmpPDB.width(6) ; // occupancy 55 - 60
06043             OutputTmpPDB.right ;
06044             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
06045             OutputTmpPDB.precision(2) ;
06046             OutputTmpPDB << dataP[i].qmGrpID ;
06047             
06048             OutputTmpPDB.width(6) ; // tempFactor/Beta 61 - 66
06049             OutputTmpPDB.right ;
06050             OutputTmpPDB.setf(std::ios::fixed,std::ios::floatfield) ;
06051             OutputTmpPDB.precision(2) ;
06052             OutputTmpPDB << dataP[i].id ;
06053             
06054             OutputTmpPDB.width(6) ; // extra information not originaly on PDB format  67 - 72
06055             OutputTmpPDB << "      " ;
06056             
06057             OutputTmpPDB.width(4) ; // segment information from NAMD, not originaly on PDB format  72 - 76
06058             OutputTmpPDB.left ;
06059             OutputTmpPDB << "QM  ";
06060             
06061             OutputTmpPDB.width(2) ; // element 77 - 78
06062             OutputTmpPDB.right ;
06063             OutputTmpPDB << element ;
06064             
06065             OutputTmpPDB.width(2) ; // charge 77 - 78
06066             OutputTmpPDB.right ;
06067             OutputTmpPDB << dataP[i].charge ;
06068             
06069             OutputTmpPDB << std::endl;
06070             
06071         }
06072         
06073         OutputTmpPDB << "END" << std::endl;
06074         
06075         OutputTmpPDB.close();
06076     }
06077     catch (...)
06078     {
06079         iout << iERROR << "Generic exception at QM write PBD." << endi ;
06080         NAMD_die("PDB write error");
06081         throw "Generic exception!" ;
06082     }
06083     return ;
06084 }
06085 
06086 #endif
06087 
06088 #include "ComputeQMMgr.def.h"
06089 

Generated on Tue Sep 25 01:17:12 2018 for NAMD by  doxygen 1.4.7