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

Generated on Sat Sep 23 01:17:13 2017 for NAMD by  doxygen 1.4.7