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