ComputeQMMgr Class Reference

List of all members.

Public Member Functions

 ComputeQMMgr ()
 ~ComputeQMMgr ()
void setCompute (ComputeQM *c)
void recvPartQM (QMCoordMsg *)
void recvFullQM (QMCoordMsg *)
void recvPntChrg (QMPntChrgMsg *)
void calcMOPAC (QMGrpCalcMsg *)
void calcORCA (QMGrpCalcMsg *)
void calcUSR (QMGrpCalcMsg *)
void storeQMRes (QMGrpResMsg *)
void procQMRes ()
void recvForce (QMForceMsg *)
SortedArray< LSSSubsDat > & get_subsArray ()

Detailed Description

Definition at line 401 of file ComputeQM.C.


Constructor & Destructor Documentation

ComputeQMMgr::ComputeQMMgr (  ) 

Definition at line 548 of file ComputeQM.C.

00548                            :
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 }

ComputeQMMgr::~ComputeQMMgr (  ) 

Definition at line 557 of file ComputeQM.C.

References close_dcd_write().

00557                             {
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 }


Member Function Documentation

void ComputeQMMgr::calcMOPAC ( QMGrpCalcMsg msg  ) 

Definition at line 2801 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, QMAtomData::charge, QMGrpCalcMsg::configLines, QMGrpCalcMsg::constants, QMGrpCalcMsg::data, DebugM, QMAtomData::dist, QMAtomData::element, endi(), QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, QMGrpCalcMsg::execPath, for(), QMGrpResMsg::force, Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), QMGrpResMsg::grpIndx, QMGrpCalcMsg::grpIndx, QMAtomData::id, iERROR(), iout, j, Vector::length(), Vector::length2(), Node::molecule, NAMD_die(), NAMD_err(), QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMGrpResMsg::numForces, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, Node::Object(), QMGrpCalcMsg::PMEEwaldCoefficient, QMGrpCalcMsg::PMEOn, QMAtomData::position, QMGrpCalcMsg::prepProc, QMGrpCalcMsg::prepProcOn, QMGrpCalcMsg::qmAtmChrgMode, QMCHRGNONE, QMPCTYPE_CLASSICAL, QMPCTYPE_IGNORE, QMGrpCalcMsg::secProc, QMGrpCalcMsg::secProcOn, QMGrpCalcMsg::switching, QMGrpCalcMsg::timestep, QMAtomData::type, Vector::unit(), QMGrpResMsg::virial, Vector::x, Vector::y, and Vector::z.

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     // For coulomb calculation
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     // We create a pointer for PME correction, which may point to
02823     // a copy of the original point charge array, with unchanged charges,
02824     // or points to the original array in case no switching or charge
02825     // scheme is being used.
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     // For each QM group, create a subdirectory where files will be palced.
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     // Opens file for coordinate and parameter input
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     // Builds the command that will be executed
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     // Builds the file name where MOPAC will place the gradient
02904     // This is also relative to the input file name.
02905     outputFileName = inputFileName ;
02906     outputFileName.append(".aux") ;
02907     
02908     if (msg->numAllPntChrgs) {
02909         // Builds the file name where we will write the point charges.
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     // writes configuration lines to the input file.
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     // write QM and dummy atom coordinates to input file and
02944     // MM electric field from MM point charges.
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             // The Electrostatic Potential is calculated according to
02959             // the QMMM Keyword documentation for MOPAC
02960             // http://openmopac.net/manual/QMMM.html
02961             pcP = msg->data + msg->numAllAtoms ;
02962             for ( size_t j=0; j < msg->numAllPntChrgs; ++j, ++pcP ) {
02963                 
02964                 // In case of QM-MM bonds, the charge of the MM1 atom is ignored
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             // We use the same Coulomb constant used in the rest of NAMD
02982             // instead of the suggested rounded 332 suggested by MOPAC.
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     // runs QM command
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     // remove coordinate file
03032 //     iret = remove(inputFileName);
03033 //     if ( iret ) { NAMD_err(0); }
03034 
03035     // remove coordinate file
03036 //     iret = remove(pntChrgFileName);
03037 //     if ( iret ) { NAMD_err(0); }
03038     
03039     // opens output file
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     // Resets the pointers.
03048     atmP = msg->data ;
03049     pcP = msg->data + msg->numAllAtoms ;
03050     
03051     // Allocates the return message.
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     // We split the data into two pointers so we can "skip" the dummy atoms
03063     // which have no representation in NAMD.
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             // We use the replace field to indicate QM atoms,
03069             // which will have charge information.
03070             resForce[i].replace = 1 ;
03071             resForce[i].id = atmP->id;
03072             atmP++;
03073         }
03074         else {
03075             // We use the replace field to indicate QM atoms,
03076             // which will have charge information.
03077             resForce[i].replace = 0 ;
03078             resForce[i].id = pcP->id;
03079             pcP++;
03080         }
03081     }
03082     
03083     // Resets the pointers.
03084     atmP = msg->data ;
03085     pcP = msg->data + msg->numAllAtoms ;
03086     
03087     // Reads the data form the output file created by the QM software.
03088     // Gradients over the QM atoms, and Charges for QM atoms will be read.
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         // We loop over the lines of the output file untill we find
03095         // the line that initiates the "atom charges" lines. Then
03096         // we read all lines and expect to get one or more charges 
03097         // per line, separated by space(s), untill we find a line that
03098         // initiates the "gradients", then once more, we expect several
03099         // numbers separated by space(s). When the "overlap matrix" 
03100         // string is found, we break the loop and stop reading the file.
03101         
03102         // if ( strstr(line,"TOTAL_ENERGY") != NULL ) {
03103         if ( strstr(line,"HEAT_OF_FORMATION") != NULL ) {
03104             
03105             char strEnergy[14], *endPtr ;
03106             
03107             //strncpy(strEnergy, line + 17, 13) ;
03108             strncpy(strEnergy, line + 28, 13) ;
03109             strEnergy[13] = '\0';
03110             
03111             // We have to convert the notation from FORTRAN double precision to
03112             // the natural, normal, reasonable and not terribly out dated format.
03113             resMsg->energyOrig = strtod(strEnergy, &endPtr);
03114             if ( *endPtr == 'D' ) {
03115                 *endPtr = 'e';
03116                 resMsg->energyOrig = strtod(strEnergy, &endPtr);
03117             }
03118             
03119             // In MOPAC, the total energy is given in EV, so we convert to Kcal/mol
03120 //             resMsg->energyOrig *= 23.060 ; // We now read Heat of Formation, which is given in Kcal/Mol
03121             
03122 //             DebugM(4,"Reading QM energy from file: " << resMsg->energyOrig << "\n");
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             // If the user asked for charges NOT to be read from MOPAC,
03135             // we skip the charge reading step.
03136             if (msg->qmAtmChrgMode == QMCHRGNONE) {
03137                 chargeFields = false;
03138                 atmIndx = msg->numAllAtoms;
03139             }
03140             else {
03141                 chargeFields = true;
03142                 atmIndx = 0;
03143             }
03144             
03145             // Now we expect the following line(s) to have atom charges
03146             continue;
03147         }
03148         
03149         if ( strstr(line,"GRADIENTS") != NULL ) {
03150             
03151             // Now that all charges have been read, checks if the 
03152             // numbers match
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             // Now we expect the following line(s) to have gradients
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             // Now that all gradients have been read, checks if the 
03171             // numbers match
03172             if (atmIndx != msg->numAllAtoms) {
03173                 NAMD_die("Error reading QM forces file. Wrong number of gradients");
03174             }
03175             
03176             gradsRead = true;
03177             
03178             // Nothing more to read from the ".aux" file
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                 // If we are reading charges from QM atoms, store them.
03194                 if (atmIndx < msg->numQMAtoms ) {
03195                     atmP[atmIndx].charge = localCharge;
03196                     resForce[atmIndx].charge = localCharge;
03197                 }
03198                 
03199                 // If we are reading charges from Dummy atoms,
03200                 // place them on the appropriate QM atoms.
03201                 if ( msg->numQMAtoms <= atmIndx &&
03202                     atmIndx < msg->numAllAtoms ) {
03203                     // We keep the calculated charge in the dummy atom
03204                     // so that we can calculate electrostatic forces later on.
03205                     atmP[atmIndx].charge = localCharge;
03206                     
03207                     // The dummy atom points to the QM atom to which it is bound.
03208                     int qmInd = atmP[atmIndx].bountToIndx ;
03209                     resForce[qmInd].charge += localCharge;
03210                 }
03211                 
03212                 strIndx += 9;
03213                 atmIndx++;
03214                 
03215                 // If we found all charges for QM and dummy atoms, break the loop
03216                 // and stop reading the line.
03217                 if (atmIndx == msg->numAllAtoms) {
03218                     chargeFields = false;
03219                     break;
03220                 }
03221                 
03222             }
03223             
03224         }
03225         
03226         int gradLength ; // Change for variable length MOPAC output
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                         // Gradients in MOPAC are written in kcal/mol/A.
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                     // If we are reading forces applied on Dummy atoms,
03251                     // place them on the appropriate QM and MM atoms to conserve energy.
03252                     
03253                     // This implementation was based on the description in 
03254                     // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
03255                     if ( msg->numQMAtoms <= atmIndx &&
03256                     atmIndx < msg->numAllAtoms ) {
03257                         // The dummy atom points to the QM atom to which it is bound.
03258                         // The QM atom and the MM atom (in a QM-MM bond) point to each other.
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                         // Unit vectors
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                 // If we found all gradients for QM atoms, break the loop
03311                 // and keep the next gradient line from being read, if there
03312                 // is one. Following gradients, if present, will refer to link
03313                 // atoms, and who cares about those?.
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     // In case charges are not to be read form the QM software output,
03328     // we load the origianl atom charges.
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             // Loops over all QM atoms (in all QM groups) comparing their global indices
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     // remove force file
03357 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
03358 //     iret = remove(outputFileName);
03359 //     if ( iret ) { NAMD_err(0); }
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     // The initial point charge index for the *force* message is the number of QM
03371     // atoms, since the dummy atoms have no representation in NAMD
03372     int pcIndx = msg->numQMAtoms;
03373     
03374     // ---- WARNING  ----
03375             // We need to apply the force felt by each QM atom due to the classical 
03376             // point charges sent to MOPAC.
03377             // MOPAC takes the MM electrostatic potential into cosideration ONLY
03378             // when performing the SCF calculation and when returning the energy
03379             // of the system, but it does not apply the potential to the final
03380             // gradient calculation, therefore, we calculate the resulting force
03381             // here.
03382     // Initialize force vector for electrostatic contribution
03383     std::vector<Force> qmElecForce ;
03384     for (size_t j=0; j<msg->numAllAtoms; ++j )
03385         qmElecForce.push_back( Force(0) );
03386     
03387     // We loop over point charges and distribute the forces applied over
03388     // virtual point charges to the MM1 and MM2 atoms (if any virtual PCs exists).
03389     for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
03390         
03391         // No force was applied to the QM region due to this charge, so it 
03392         // does not receive any back from the QM region. It must be an MM1 atom
03393         // from a QM-MM bond.
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             // We accumulate the total force felt by a point charge
03414             // due to the QM charge distribution. This total force
03415             // will be applied to the point charge if it is a real one,
03416             // or will be distirbuted over MM1 and MM2 point charges, it 
03417             // this is a virtual point charge.
03418             totalForce += force*rVec.unit();
03419             
03420             // Accumulate force felt by a QM atom due to point charges.
03421             qmElecForce[j] += -1*force*rVec.unit();
03422         }
03423         
03424         if (pcP[i].type == QMPCTYPE_CLASSICAL) {
03425             // Checking pcP was not a QM atom in another region
03426             // If so the interaction is accounted in the other region
03427             if (qmAtomGroup[pcP[i].id] == 0) {
03428             // If we already ignored MM1 charges, we take all other 
03429             // non-virtual charges and apply forces directly to them.
03430                 resForce[pcIndx].force += totalForce;
03431             }
03432         }
03433         else {
03434             // If we are handling virtual PC, we distribute the force over
03435             // MM1 and MM2.
03436             
03437             // Virtual PC are bound to MM2.
03438             int mm2Indx = pcP[i].bountToIndx;
03439             // MM2 charges are bound to MM1.
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     // We now apply the accumulated electrostatic force contribution
03454     // to QM atoms.
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             // In case the QM atom is a Link atom, we redistribute 
03462             // its force as bove.
03463             
03464             // The dummy atom points to the QM atom to which it is bound.
03465             // The QM atom and the MM atom (in a QM-MM bond) point to each other.
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             // Unit vectors
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     // Adjusts forces from PME, canceling contributions from the QM and 
03509     // direct Coulomb forces calculated here.
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 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
03539                 BigReal recip_energy = (1-tmp_b)/r;
03540                 
03541 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
03542                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
03543                 
03544                 // Final force and energy correction for this pair of atoms.
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                 // The force is *subtracted* from the total force acting on
03550                 // both atoms. The sign on fixForce corrects the orientation
03551                 // of the subtracted force.
03552 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
03553 //                     << std::endl);
03554 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
03555 //                     << std::endl);
03556 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
03557 //                 DebugM(4,"Energy correction: " << energy << "\n");
03558                 resForce[i].force -= fixForce ;
03559                 resForce[j].force -= -1*fixForce ;
03560                 
03561                 // The energy is *subtracted* from the total energy calculated here.
03562                 resMsg->energyCorr -= energy;
03563             }
03564         }
03565         
03566 //         DebugM(4,"Corrected energy QM-QM interactions: " << resMsg->energyCorr << "\n");
03567         
03568         pcIndx = msg->numQMAtoms;
03569         // We only loop over point charges from real atoms, ignoring the ones 
03570         // created to handle QM-MM bonds.
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 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
03594                 BigReal recip_energy = (1-tmp_b)/r;
03595                 
03596 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
03597                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
03598                 
03599                 // Final force and energy correction for this pair of atoms.
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                 // The energy is *subtracted* from the total energy calculated here.
03605                 resMsg->energyCorr -= energy;
03606             }
03607             
03608             // The force is *subtracted* from the total force acting on
03609             // the point charge..
03610 //             DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
03611 //             << std::endl);
03612 //             DebugM(4,"Force correction: " << fixForce << std::endl);
03613             resForce[pcIndx].force -= kq_i*fixForce ;
03614         }
03615         
03616     }
03617     
03618     // Calculates the virial contribution form the forces on QM atoms and 
03619     // point charges calculated here.
03620     for (size_t i=0; i < msg->numQMAtoms; i++) {
03621         // virial used by NAMD is -'ve of normal convention, so reverse it!
03622         // virial[i][j] in file should be sum of -1 * f_i * r_j
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; // Index in the real PC force array.
03629     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
03630         // virial used by NAMD is -'ve of normal convention, so reverse it!
03631         // virial[i][j] in file should be sum of -1 * f_i * r_j
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     // Send message to rank zero with results.
03640     QMProxy[0].recvQMRes(resMsg);
03641     
03642     if (msg->switching && msg->PMEOn)
03643         delete [] pcPpme;
03644     
03645     delete msg;
03646     return ;
03647 }

void ComputeQMMgr::calcORCA ( QMGrpCalcMsg msg  ) 

Definition at line 3649 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, QMAtomData::charge, QMGrpCalcMsg::charge, QMGrpCalcMsg::configLines, QMGrpCalcMsg::constants, QMGrpCalcMsg::data, DebugM, QMAtomData::dist, endi(), QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, QMGrpCalcMsg::execPath, for(), QMGrpResMsg::force, Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), QMGrpResMsg::grpIndx, QMGrpCalcMsg::grpIndx, QMAtomData::id, iERROR(), iout, j, Vector::length(), Node::molecule, QMGrpCalcMsg::multiplicity, NAMD_die(), NAMD_err(), QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMGrpResMsg::numForces, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, Node::Object(), QMGrpCalcMsg::PMEEwaldCoefficient, QMGrpCalcMsg::PMEOn, QMAtomData::position, QMGrpCalcMsg::prepProc, QMGrpCalcMsg::prepProcOn, QMGrpCalcMsg::qmAtmChrgMode, QMCHRGCHELPG, QMCHRGMULLIKEN, QMCHRGNONE, QMPCTYPE_CLASSICAL, QMPCTYPE_IGNORE, QMGrpCalcMsg::secProc, QMGrpCalcMsg::secProcOn, QMGrpCalcMsg::switching, QMGrpCalcMsg::timestep, QMAtomData::type, QMGrpResMsg::virial, Vector::x, Vector::y, and Vector::z.

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     // For coulomb calculation
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     // We create a pointer for PME correction, which may point to
03671     // a copy of the original point charge array, with unchanged charges,
03672     // or points to the original array in case no switching or charge
03673     // scheme is being used.
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     // For each QM group, create a subdirectory where files will be palced.
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     // Opens file for coordinate and parameter input
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     // Builds the command that will be executed
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     // Build the redirect file name we need for the screen output.
03751     // That's the only place where we can get partial charges for QM atoms.
03752     tmpRedirectFileName = inputFileName ;
03753     tmpRedirectFileName.append(".TmpOut") ;
03754     
03755     qmCommand.append(" > ") ;
03756     qmCommand.append(tmpRedirectFileName) ;
03757     
03758     // Builds the file name where orca will place the gradient
03759     // This will be relative to the input file
03760     outputFileName = inputFileName ;
03761     outputFileName.append(".engrad") ;
03762     
03763     // Builds the file name where orca will place gradients acting on 
03764     // point charges
03765     pcGradFilename = inputFileName ;
03766     pcGradFilename.append(".pcgrad") ;
03767     
03768     if (msg->numAllPntChrgs) {
03769         // Builds the file name where we will write the point charges.
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     // writes configuration lines to the input file.
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     // write QM and dummy atom coordinates to input file.
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         // Write point charges to file.
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         // runs QM command
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     // remove coordinate file
03899 //     iret = remove(inputFileName);
03900 //     if ( iret ) { NAMD_err(0); }
03901 
03902     // remove coordinate file
03903 //     iret = remove(pntChrgFileName);
03904 //     if ( iret ) { NAMD_err(0); }
03905     
03906     // opens output file
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     // Resets the pointers.
03915     atmP = msg->data ;
03916     pcP = msg->data + msg->numAllAtoms ;
03917     
03918     // Allocates the return message.
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     // We split the data into two pointers so we can "skip" the dummy atoms
03930     // which have no representation in NAMD.
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             // We use the replace field to indicate QM atoms,
03936             // which will have charge information.
03937             resForce[i].replace = 1 ;
03938             resForce[i].id = atmP->id;
03939             atmP++;
03940         }
03941         else {
03942             // We use the replace field to indicate QM atoms,
03943             // which will have charge information.
03944             resForce[i].replace = 0 ;
03945             resForce[i].id = pcP->id;
03946             pcP++;
03947         }
03948     }
03949     
03950     // Resets the pointers.
03951     atmP = msg->data ;
03952     pcP = msg->data + msg->numAllAtoms ;
03953     
03954     size_t atmIndx = 0, gradCount = 0;
03955     int numAtomsInFile = 0;
03956     
03957     // Reads the data form the output file created by the QM software.
03958     // Gradients over the QM atoms, and Charges for QM atoms will be read.
03959     
03960     // skip lines before number of atoms
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     // skip lines before energy
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 //     iout << "Energy in step (Hartree): " << resMsg->energyOrig << "\n" <<  endi;
03987     // All energies are given in Eh (Hartree)
03988     // NAMD needs energies in kcal/mol
03989     // The conversion factor is 627.509469
03990     resMsg->energyOrig *= 627.509469;
03991 //     iout << "Energy in step (Kcal/mol): " << resMsg->energyOrig << "\n" <<  endi;
03992     
03993     resMsg->energyCorr = resMsg->energyOrig;
03994     
03995     // skip lines before gradient
03996     for (int i = 0; i < 3; i++) {
03997         fgets(line, lineLen, outputFile) ;
03998     }
03999     
04000     // Break the loop when we find all gradients for QM atoms, 
04001     // and keep the next gradient lines from being read, if there
04002     // are more. Following gradients, if present, will refer to link
04003     // atoms.
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             // All gradients are given in Eh/a0 (Hartree over Bohr radius)
04014             // NAMD needs forces in kcal/mol/angstrons
04015             // The conversion factor is 627.509469/0.529177 = 1185.82151
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             // If we are reading forces applied on Dummy atoms,
04023             // place them on the appropriate QM and MM atoms to conserve energy.
04024             
04025             // This implementation was based on the description in 
04026             // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
04027             if ( msg->numQMAtoms <= atmIndx &&
04028             atmIndx < msg->numAllAtoms ) {
04029                 // The dummy atom points to the QM atom to which it is bound.
04030                 // The QM atom and the MM atom (in a QM-MM bond) point to each other.
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                 // Unit vectors
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     // In case charges are not to be read form the QM software output,
04082     // we load the origianl atom charges.
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             // Loops over all QM atoms (in all QM groups) comparing their global indices
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         // opens redirected output file
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             // We loop over the lines of the output file untill we find
04122             // the line that initiates the charges lines. Then
04123             // we read all lines and expect to get one charge
04124             // per line, untill we find a line that sums all charges.
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                     // Now we expect the following line to have a series of dashes
04136                     // and the folowing lines to have atom charges (among other info)
04137                     continue;
04138                 }
04139                 
04140                 if ( strstr(line,"Sum of atomic charges") != NULL ) {
04141                     
04142                     // Now that all charges have been read, grabs the total charge
04143                     // just for fun... and checking purposes.
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                     // Nothing more to read from the output file
04153                     break;
04154                 }
04155                 
04156                 // Line state 1 means we have to skip ONLY one line.
04157                 if (lineState == 1) {
04158                     lineState = 2;
04159                     continue;
04160                 }
04161                 
04162                 // Line state 2 means we have to read the line and grab the info.
04163                 if (lineState == 2) {
04164                     
04165                     strncpy(result, line+8,12) ;
04166                     result[12] = '\0';
04167                     
04168                     Real localCharge = atof(result);
04169                     
04170                     // If we are reading charges from QM atoms, store them.
04171                     if (atmIndx < msg->numQMAtoms ) {
04172                         atmP[atmIndx].charge = localCharge;
04173                         resForce[atmIndx].charge = localCharge;
04174                     }
04175                     
04176                     // If we are reading charges from Dummy atoms,
04177                     // place the on the appropriate QM atom.
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                     // If we found all charges for QM atoms, change the lineState
04188                     // untill we reach the "total charge" line.
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                     // Now we expect the following line to have a series of dashes
04206                     // and the folowing lines to have atom charges (among other info)
04207                     continue;
04208                 }
04209                 
04210                 if ( strstr(line,"Total charge") != NULL ) {
04211                     
04212                     // Now that all charges have been read, grabs the total charge
04213                     // just for fun... and checking purposes.
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                     // Nothing more to read from the output file
04223                     break;
04224                 }
04225                 
04226                 // Line state 1 means we have to skip ONLY one line.
04227                 if (lineState == 1) {
04228                     lineState = 2;
04229                     continue;
04230                 }
04231                 
04232                 // Line state 2 means we have to read the line and grab the info.
04233                 if (lineState == 2) {
04234                     
04235                     strncpy(result, line+12,15) ;
04236                     result[15] = '\0';
04237                     
04238                     Real localCharge = atof(result);
04239                     
04240                     // If we are reading charges from QM atoms, store them.
04241                     if (atmIndx < msg->numQMAtoms ) {
04242                         atmP[atmIndx].charge = localCharge;
04243                         resForce[atmIndx].charge = localCharge;
04244                     }
04245                     
04246                     // If we are reading charges from Dummy atoms,
04247                     // place the on the appropriate QM atom.
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                     // If we found all charges for QM atoms, we ignore the following line
04258                     // untill we reach the "total charge" line.
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     // remove force file
04274 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
04275 //     iret = remove(outputFileName);
04276 //     if ( iret ) { NAMD_err(0); }
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     // The initial point charge index for the force message is the number of QM
04285     // atoms, since the dummy atoms have no representation in NAMD
04286     int pcIndx = msg->numQMAtoms;
04287     
04288     // If there are no point charges, orca will not create a "pcgrad" file.
04289     if (msg->numAllPntChrgs > 0) {
04290     
04291         // opens redirected output file
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         // We loop over point charges, reading the total electrostatic force 
04309         // applied on them by the QM region.
04310         // We redistribute the forces applied over virtual point
04311         // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
04312         for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
04313             
04314             Force totalForce(0);
04315             
04316             // No force was applied to the QM region due to this charge, since it
04317             // was skipped when writing down point charges to the QM software, so it
04318             // does not receive any back from the QM region. It must be an MM1 atom
04319             // from a QM-MM bond.
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             // All gradients are given in Eh/a0 (Hartree over Bohr radius)
04330             // NAMD needs forces in kcal/mol/angstrons
04331             // The conversion factor is 627.509469/0.529177 = 1185.82151
04332             totalForce *= -1185.82151;
04333             
04334             if (pcP[i].type == QMPCTYPE_CLASSICAL) {
04335                 // Checking pcP was not a QM atom in another region
04336                 // If so the interaction is accounted in the other region
04337                 if (qmAtomGroup[pcP[i].id] == 0) {
04338                 // If we already ignored MM1 charges, we take all other 
04339                 // non-virtual charges and apply forces directly to them.
04340                     resForce[pcIndx].force += totalForce;
04341                 }
04342             }
04343             else {
04344                 // If we are handling virtual PC, we distribute the force over
04345                 // MM1 and MM2.
04346                 
04347                 Force mm1Force(0), mm2Force(0);
04348                 
04349                 // Virtual PC are bound to MM2.
04350                 int mm2Indx = pcP[i].bountToIndx;
04351                 // MM2 charges are bound to MM1.
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     // Adjusts forces from PME, canceling contributions from the QM and 
04371     // direct Coulomb forces calculated here.
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 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
04399                 BigReal recip_energy = (1-tmp_b)/r;
04400                 
04401 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
04402                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
04403                 
04404                 const BigReal kq_i = p_i_charge * constants;
04405                 
04406                 // Final force and energy correction for this pair of atoms.
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                 // The force is *subtracted* from the total force acting on
04412                 // both atoms. The sign on fixForce corrects the orientation
04413                 // of the subtracted force.
04414 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
04415 //                     << std::endl);
04416 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
04417 //                     << std::endl);
04418 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
04419                 resForce[i].force -= fixForce ;
04420                 resForce[j].force -= -1*fixForce;
04421                 
04422                 // The energy is *subtracted* from the total energy calculated here.
04423                 resMsg->energyCorr -= energy;
04424             }
04425         }
04426         
04427         pcIndx = msg->numQMAtoms;
04428         // We only loop over point charges from real atoms, ignoring the ones 
04429         // created to handle QM-MM bonds.
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 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
04453                 BigReal recip_energy = (1-tmp_b)/r;
04454                 
04455 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
04456                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
04457                 
04458                 // Final force and energy correction for this pair of atoms.
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                 // The energy is *subtracted* from the total energy calculated here.
04464                 resMsg->energyCorr -= energy;
04465                 
04466             }
04467             
04468             // The force is *subtracted* from the total force acting on
04469                 // the point charge.
04470 //                 DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
04471 //                     << std::endl);
04472 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
04473             resForce[pcIndx].force -= kq_i*fixForce ;
04474         }
04475         
04476     }
04477     
04478     DebugM(1,"Determining virial...\n");
04479     
04480     // Calculates the virial contribution form the forces on QM atoms and 
04481     // point charges calculated here.
04482     for (size_t i=0; i < msg->numQMAtoms; i++) {
04483         // virial used by NAMD is -'ve of normal convention, so reverse it!
04484         // virial[i][j] in file should be sum of -1 * f_i * r_j
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; // Index in the real PC force array.
04491     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
04492         // virial used by NAMD is -'ve of normal convention, so reverse it!
04493         // virial[i][j] in file should be sum of -1 * f_i * r_j
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     // Send message to rank zero with results.
04502     QMProxy[0].recvQMRes(resMsg);
04503     
04504     if (msg->switching && msg->PMEOn)
04505         delete [] pcPpme;
04506     
04507     delete msg;
04508     return ;
04509 }

void ComputeQMMgr::calcUSR ( QMGrpCalcMsg msg  ) 

Definition at line 4511 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, QMAtomData::charge, QMGrpCalcMsg::constants, QMGrpCalcMsg::data, DebugM, QMAtomData::dist, endi(), QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, QMGrpCalcMsg::execPath, for(), QMGrpResMsg::force, Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), QMGrpResMsg::grpIndx, QMGrpCalcMsg::grpIndx, QMAtomData::id, iERROR(), iout, j, Vector::length(), Vector::length2(), Node::molecule, NAMD_die(), NAMD_err(), QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMGrpResMsg::numForces, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, Node::Object(), QMGrpCalcMsg::PMEEwaldCoefficient, QMGrpCalcMsg::PMEOn, QMAtomData::position, QMGrpCalcMsg::prepProc, QMGrpCalcMsg::prepProcOn, QMGrpCalcMsg::qmAtmChrgMode, QMCHRGNONE, QMPCTYPE_CLASSICAL, QMPCTYPE_IGNORE, QMGrpCalcMsg::secProc, QMGrpCalcMsg::secProcOn, QMGrpCalcMsg::switching, QMGrpCalcMsg::timestep, QMAtomData::type, Vector::unit(), QMGrpResMsg::virial, Vector::x, Vector::y, and Vector::z.

04511                                             {
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     // For coulomb calculation
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     // We create a pointer for PME correction, which may point to
04526     // a copy of the original point charge array, with unchanged charges,
04527     // or points to the original array in case no switching or charge
04528     // scheme is being used.
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     // For each QM group, create a subdirectory where files will be palced.
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     // Opens file for coordinate and parameter input
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     // Builds the command that will be executed
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     // Builds the file name where orca will place the gradient
04606     // This will be relative to the input file
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     // write QM and dummy atom coordinates to input file.
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     // Write point charges to file.
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         // runs QM command
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     // remove coordinate file
04695 //     iret = remove(inputFileName);
04696 //     if ( iret ) { NAMD_err(0); }
04697 
04698     // remove coordinate file
04699 //     iret = remove(pntChrgFileName);
04700 //     if ( iret ) { NAMD_err(0); }
04701     
04702     // opens output file
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     // Resets the pointers.
04711     atmP = msg->data ;
04712     pcP = msg->data + msg->numAllAtoms ;
04713     
04714     // Allocates the return message.
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     // We split the data into two pointers so we can "skip" the dummy atoms
04726     // which have no representation in NAMD.
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             // We use the replace field to indicate QM atoms,
04732             // which will have charge information.
04733             resForce[i].replace = 1 ;
04734             resForce[i].id = atmP->id;
04735             atmP++;
04736         }
04737         else {
04738             // We use the replace field to indicate QM atoms,
04739             // which will have charge information.
04740             resForce[i].replace = 0 ;
04741             resForce[i].id = pcP->id;
04742             pcP++;
04743         }
04744     }
04745     
04746     // Resets the pointers.
04747     atmP = msg->data ;
04748     pcP = msg->data + msg->numAllAtoms ;
04749     
04750     // Reads the data form the output file created by the QM software.
04751     // Gradients over the QM atoms, and Charges for QM atoms will be read.
04752     
04753     // Number of point charges for which we will receive forces.
04754     int usrPCnum = 0;
04755     const size_t lineLen = 256;
04756     char *line = new char[lineLen];
04757     
04758     fgets(line, lineLen, outputFile);
04759     
04760 //     iret = fscanf(outputFile,"%lf %d\n", &resMsg->energyOrig, &usrPCnum);
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         // If we are reading charges and forces on QM atoms, store
04788         // them directly.
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         // If we are reading charges and forces applied on Dummy atoms,
04800         // place them on the appropriate QM and MM atoms to conserve energy.
04801         
04802         // This implementation was based on the description in 
04803         // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
04804         if ( msg->numQMAtoms <= atmIndx &&
04805         atmIndx < msg->numAllAtoms ) {
04806             
04807             // We keep the calculated charge in the dummy atom
04808             // so that we can calculate electrostatic forces later on.
04809             atmP[atmIndx].charge = localCharge;
04810             
04811             // If we are reading charges from Dummy atoms,
04812             // place them on the appropriate QM atom.
04813             // The dummy atom points to the QM atom to which it is bound.
04814             int qmInd = atmP[atmIndx].bountToIndx ;
04815             resForce[qmInd].charge += localCharge;
04816             
04817             // The dummy atom points to the QM atom to which it is bound.
04818             // The QM atom and the MM atom (in a QM-MM bond) point to each other.
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             // Unit vectors
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     // The initial point charge index for the force message is the number of QM
04862     // atoms, since the dummy atoms have no representation in NAMD
04863     int pcIndx = msg->numQMAtoms;
04864     
04865     if (usrPCnum > 0) {
04866         // We loop over point charges, reading the total electrostatic force 
04867         // applied on them by the QM region.
04868         // We redistribute the forces applied over virtual point
04869         // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
04870         for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
04871             
04872             Force totalForce(0);
04873             
04874             // No force was applied to the QM region due to this charge, since it
04875             // was skipped when writing down point charges to the QM software, so it
04876             // does not receive any back from the QM region. It must be an MM1 atom
04877             // from a QM-MM bond.
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                 // Checking pcP was not a QM atom in another region
04889                 // If so the interaction is accounted in the other region
04890                 if (qmAtomGroup[pcP[i].id] == 0) {
04891                 // If we already ignored MM1 charges, we take all other 
04892                 // non-virtual charges and apply forces directly to them.
04893                     resForce[pcIndx].force += totalForce;
04894                 }
04895             }
04896             else {
04897                 // If we are handling virtual PC, we distribute the force over
04898                 // MM1 and MM2.
04899                 
04900                 Force mm1Force(0), mm2Force(0);
04901                 
04902                 // Virtual PC are bound to MM2.
04903                 int mm2Indx = pcP[i].bountToIndx;
04904                 // MM2 charges are bound to MM1.
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     // In case charges are not to be read form the QM software output,
04924     // we load the origianl atom charges.
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             // Loops over all QM atoms (in all QM groups) comparing their global indices
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     // remove force file
04953 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
04954 //     iret = remove(outputFileName);
04955 //     if ( iret ) { NAMD_err(0); }
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         // We only loop over point charges from real atoms, ignoring the ones 
04964         // created to handle QM-MM bonds.
04965         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
04966             
04967             // No force was applied to the QM region due to this charge, so it 
04968             // does not receive any back from the QM region. It must be an MM1 atom
04969             // from a QM-MM bond.
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                 // We accumulate the total force felt by a point charge
04990                 // due to the QM charge distribution. This total force
04991                 // will be applied to the point charge if it is a real one,
04992                 // or will be distirbuted over MM1 and MM2 point charges, it 
04993                 // this is a virtual point charge.
04994                 totalForce += force*rVec.unit();
04995             }
04996             
04997             if (pcP[i].type == QMPCTYPE_CLASSICAL) {
04998                 // Checking pcP was not a QM atom in another region
04999                 // If so the interaction is accounted in the other region
05000                 if (qmAtomGroup[pcP[i].id] == 0) {
05001                 // If we already ignored MM1 charges, we take all other 
05002                 // non-virtual charges and apply forces directly to them.
05003                     resForce[pcIndx].force += totalForce;
05004                 }
05005             }
05006             else {
05007                 // If we are handling virtual PC, we distribute the force over
05008                 // MM1 and MM2.
05009                 
05010                 Force mm1Force(0), mm2Force(0);
05011                 
05012                 // Virtual PC are bound to MM2.
05013                 int mm2Indx = pcP[i].bountToIndx;
05014                 // MM2 charges are bound to MM1.
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     // Adjusts forces from PME, canceling contributions from the QM and 
05030     // direct Coulomb forces calculated here.
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 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
05058                 BigReal recip_energy = (1-tmp_b)/r;
05059                 
05060 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
05061                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
05062                 
05063                 const BigReal kq_i = p_i_charge * constants;
05064                 
05065                 // Final force and energy correction for this pair of atoms.
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                 // The force is *subtracted* from the total force acting on
05071                 // both atoms. The sign on fixForce corrects the orientation
05072                 // of the subtracted force.
05073 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
05074 //                     << std::endl);
05075 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
05076 //                     << std::endl);
05077 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
05078                 resForce[i].force -= fixForce ;
05079                 resForce[j].force -= -1*fixForce;
05080                 
05081                 // The energy is *subtracted* from the total energy calculated here.
05082                 resMsg->energyCorr -= energy;
05083             }
05084         }
05085         
05086         pcIndx = msg->numQMAtoms;
05087         // We only loop over point charges from real atoms, ignoring the ones 
05088         // created to handle QM-MM bonds.
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 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
05112                 BigReal recip_energy = (1-tmp_b)/r;
05113                 
05114 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
05115                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
05116                 
05117                 // Final force and energy correction for this pair of atoms.
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                 // The energy is *subtracted* from the total energy calculated here.
05123                 resMsg->energyCorr -= energy;
05124                 
05125             }
05126             
05127             // The force is *subtracted* from the total force acting on
05128                 // the point charge.
05129 //                 DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
05130 //                     << std::endl);
05131 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
05132             resForce[pcIndx].force -= kq_i*fixForce ;
05133         }
05134         
05135     }
05136     
05137     DebugM(1,"Determining virial...\n");
05138     
05139     // Calculates the virial contribution form the forces on QM atoms and 
05140     // point charges calculated here.
05141     for (size_t i=0; i < msg->numQMAtoms; i++) {
05142         // virial used by NAMD is -'ve of normal convention, so reverse it!
05143         // virial[i][j] in file should be sum of -1 * f_i * r_j
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; // Index in the real PC force array.
05150     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
05151         // virial used by NAMD is -'ve of normal convention, so reverse it!
05152         // virial[i][j] in file should be sum of -1 * f_i * r_j
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     // Send message to rank zero with results.
05161     QMProxy[0].recvQMRes(resMsg);
05162     
05163     if (msg->switching && msg->PMEOn)
05164         delete [] pcPpme;
05165     
05166     delete msg;
05167     return ;
05168 }

SortedArray<LSSSubsDat>& ComputeQMMgr::get_subsArray (  )  [inline]

Definition at line 428 of file ComputeQM.C.

Referenced by lssSubs().

00428 { return subsArray; } ;

void ComputeQMMgr::procQMRes (  ) 

Definition at line 2485 of file ComputeQM.C.

References QMForce::charge, QMPntChrgMsg::coord, QMCoordMsg::coord, DebugM, endi(), QMForceMsg::energy, QMForceMsg::force, Molecule::get_numQMAtoms(), ComputeQMPntChrg::id, ComputeQMAtom::id, iINFO(), SortedArray< Elem >::insert(), iout, j, QMForceMsg::lssDat, QMPntChrgMsg::numAtoms, QMCoordMsg::numAtoms, QMForceMsg::numForces, QMForceMsg::numLssDat, SimParameters::PMEOn, QMForceMsg::PMEOn, ComputeQMAtom::position, SimParameters::qmOutFreq, SimParameters::qmPosOutFreq, ResizeArray< Elem >::size(), SortedArray< Elem >::sort(), QMPntChrgMsg::sourceNode, QMCoordMsg::sourceNode, QMForceMsg::virial, write_dcdstep(), Vector::x, Vector::y, and Vector::z.

02485                              {
02486     
02487     // Writes a DCD file with the charges of all QM atoms at a frequency 
02488     // defined by the user in qmOutFreq.
02489     if ( simParams->qmOutFreq > 0 && 
02490          timeStep % simParams->qmOutFreq == 0 ) {
02491         
02492         iout << iINFO << "Writing QM charge output at step " 
02493         << timeStep <<  "\n" << endi;
02494     
02495         Real *x = outputData, 
02496              *y = outputData + molPtr->get_numQMAtoms(), 
02497              *z = outputData + 2*molPtr->get_numQMAtoms();
02498         
02499         for (int qmIt=0; qmIt<numQMAtms; qmIt++){
02500             x[qmIt] = qmCoord[qmIt].id;
02501             y[qmIt] = force[qmCoord[qmIt].id].charge ;
02502             z[qmIt] = 0;
02503         }
02504         
02505         write_dcdstep(dcdOutFile, numQMAtms, x, y, z, 0) ;
02506     }
02507     
02508     // Writes a DCD file with the positions of all QM atoms at a frequency 
02509     // defined by the user in qmPosOutFreq.
02510     if ( simParams->qmPosOutFreq > 0 && 
02511          timeStep % simParams->qmPosOutFreq == 0 ) {
02512         
02513         iout << iINFO << "Writing QM position output at step " 
02514         << timeStep <<  "\n" << endi;
02515         
02516         SortedArray<idIndxStr> idIndx;
02517         
02518         for(int i=0; i<numQMAtms;i++) {
02519             idIndx.insert( idIndxStr(qmCoord[i].id, i) );
02520         }
02521         idIndx.sort();
02522         
02523         Real *x = outputData, 
02524              *y = outputData + molPtr->get_numQMAtoms(), 
02525              *z = outputData + 2*molPtr->get_numQMAtoms();
02526         
02527         for (int qmIt=0; qmIt<numQMAtms; qmIt++){
02528             x[qmIt] = qmCoord[idIndx[qmIt].indx].position.x;
02529             y[qmIt] = qmCoord[idIndx[qmIt].indx].position.y;
02530             z[qmIt] = qmCoord[idIndx[qmIt].indx].position.z;
02531         }
02532         
02533         write_dcdstep(dcdPosOutFile, numQMAtms, x, y, z, 0) ;
02534     }
02535     
02536     // distribute forces
02537     DebugM(4,"Distributing QM forces for all ranks.\n");
02538     for ( int j=0; j < numSources; ++j ) {
02539         
02540         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             // Since we receive two different messages from nodes, there is no 
02565             // guarantee the two sets of messages will come in the same order.
02566             // Therefore, we match the messages by comaring their sourceNodes.
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             // I want to count number of different atom ids
02586             // which is the size of force array. 
02587             // Avoids double counting of forces
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(); // Fixes same patch different QM regions double counting
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         // This iterator is used in BOTH loops!
02605         int forceIter = 0;
02606         auxset0.clear(); // Need to keep track of forces by id that have been copied to fmsg
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); // This should add each qm atom once
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                 // (QM atoms that are PC or repeated PC) are not copied to fmsg
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 }

void ComputeQMMgr::recvForce ( QMForceMsg fmsg  ) 

Definition at line 2662 of file ComputeQM.C.

References SortedArray< Elem >::add(), QMForceMsg::lssDat, QMForceMsg::numLssDat, and ComputeQM::saveResults().

02662                                              {
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 }

void ComputeQMMgr::recvFullQM ( QMCoordMsg qmFullMsg  ) 

Definition at line 1281 of file ComputeQM.C.

References ResizeArray< Elem >::clear(), ComputeQM::processFullQM(), and ResizeArray< Elem >::size().

01281                                                    {
01282     
01283     if (subsArray.size() > 0)
01284         subsArray.clear();
01285     
01286     QMCompute->processFullQM(qmFullMsg);
01287 }

void ComputeQMMgr::recvPartQM ( QMCoordMsg msg  ) 

Definition at line 818 of file ComputeQM.C.

References ResizeArray< Elem >::add(), QMForce::charge, ComputeQMAtom::charge, ComputeQMPntChrg::charge, QMPntChrgMsg::coord, QMCoordMsg::coord, DCD_FILEEXISTS, DebugM, ComputeQMPntChrg::dist, SimParameters::dt, endi(), SimParameters::firstTimestep, QMForce::force, Molecule::get_cSMDcoffs(), Molecule::get_cSMDindex(), Molecule::get_cSMDindxLen(), Molecule::get_cSMDKs(), Molecule::get_cSMDnumInst(), Molecule::get_cSMDpairs(), Molecule::get_cSMDVels(), Molecule::get_numQMAtoms(), Molecule::get_qmcSMD(), Molecule::get_qmGrpChrg(), Molecule::get_qmGrpID(), Molecule::get_qmGrpMult(), Molecule::get_qmLSSFreq(), Molecule::get_qmLSSIdxs(), Molecule::get_qmLSSMass(), Molecule::get_qmLSSRefIDs(), Molecule::get_qmLSSRefSize(), Molecule::get_qmLSSResSize(), Molecule::get_qmLSSSize(), Molecule::get_qmMeNumBonds(), Molecule::get_qmNumGrps(), Molecule::get_qmPCFreq(), Molecule::get_qmReplaceAll(), QMForce::homeIndx, ComputeQMAtom::homeIndx, ComputeQMPntChrg::homeIndx, QMForce::id, ComputeQMAtom::id, ComputeQMPntChrg::id, iERROR(), iINFO(), iout, iWARN(), j, ComputeQMPntChrg::mass, Node::molecule, NAMD_bug(), NAMD_die(), NAMD_err(), QMPntChrgMsg::numAtoms, QMCoordMsg::numAtoms, Molecule::numAtoms, QMCoordMsg::numPCIndxs, Node::Object(), PatchMap::Object(), open_dcd_write(), SimParameters::outputFilename, QMCoordMsg::pcIndxs, PCMODECUSTOMSEL, PCMODEUPDATEPOS, PCMODEUPDATESEL, QMCoordMsg::pcSelMode, ComputeQMAtom::position, ComputeQMPntChrg::position, SimParameters::qmBondValType, SimParameters::qmCustomPCSel, ComputeQMAtom::qmGrpID, ComputeQMPntChrg::qmGrpID, SimParameters::qmLSSOn, SimParameters::qmNoPC, SimParameters::qmOutFreq, SimParameters::qmPosOutFreq, SimParameters::qmSimsPerNode, QMForce::replace, Node::simParameters, ResizeArray< Elem >::size(), QMPntChrgMsg::sourceNode, QMCoordMsg::sourceNode, TIMEFACTOR, QMCoordMsg::timestep, ComputeQMPntChrg::vdwType, and write_dcdheader().

00819 {
00820     // In the first (ever) step of the simulation, we allocate arrays that
00821     // will be used throughout the simulation, and get some important numbers.
00822     if ( ! numSources ) {
00823         DebugM(4,"Initializing ComputeQMMgr variables." << std::endl);
00824         numSources = (PatchMap::Object())->numNodesWithPatches();
00825         
00826         DebugM(4,"There are " << numSources << " nodes with patches." << std::endl);
00827         qmCoordMsgs = new QMCoordMsg*[numSources];
00828         for ( int i=0; i<numSources; ++i ) { 
00829             qmCoordMsgs[i] = 0;
00830         }
00831         
00832         // Prepares the allocation for the recvPntChrg function.
00833         
00834         DebugM(4,"Getting data from molecule and simParameters." << std::endl);
00835         
00836         molPtr = Node::Object()->molecule;
00837         simParams = Node::Object()->simParameters;
00838         
00839         numAtoms = molPtr->numAtoms;
00840         force = new QMForce[numAtoms];
00841         
00842         numQMAtms = molPtr->get_numQMAtoms();
00843         qmAtmIter = 0;
00844         
00845         noPC = simParams->qmNoPC ;
00846         meNumMMIndx = molPtr->get_qmMeNumBonds();
00847         if (noPC && meNumMMIndx == 0) {
00848             pntChrgCoordMsgs = NULL;
00849         }
00850         else {
00851             pntChrgCoordMsgs = new QMPntChrgMsg*[numSources];
00852             for ( int i=0; i<numSources; ++i ) { 
00853                 pntChrgCoordMsgs[i] = 0;
00854             }
00855         }
00856         
00857         qmPCFreq = molPtr->get_qmPCFreq();
00858         resendPCList = false;
00859         
00860         grpID = molPtr->get_qmGrpID() ;
00861         bondValType = simParams->qmBondValType;
00862         
00863         numQMGrps = molPtr->get_qmNumGrps();
00864         
00865         grpChrg = molPtr->get_qmGrpChrg() ;
00866         
00867         grpMult = molPtr->get_qmGrpMult() ;
00868         
00869         qmLSSOn = simParams->qmLSSOn ;
00870         if (qmLSSOn) {
00871             qmLSSFreq = molPtr->get_qmLSSFreq() ;
00872             qmLSSSize = molPtr->get_qmLSSSize() ;
00873             qmLSSIdxs = molPtr->get_qmLSSIdxs() ;
00874             qmLSSMass = molPtr->get_qmLSSMass() ;
00875             qmLSSResSize = molPtr->get_qmLSSResSize() ;
00876             qmLSSRefIDs = molPtr->get_qmLSSRefIDs() ;
00877             qmLSSRefSize = molPtr->get_qmLSSRefSize() ;
00878             
00879             lssPrepare();
00880         }
00881         
00882         numArrivedQMMsg = 0 ;
00883         numArrivedPntChrgMsg = 0 ;
00884         
00885         qmCoord = new ComputeQMAtom[numQMAtms];
00886         
00887         replaceForces = 0;
00888         if (molPtr->get_qmReplaceAll()) {
00889             replaceForces = 1;
00890         }
00891         
00892         cSMDon = molPtr->get_qmcSMD() ;
00893         if (cSMDon) {
00894                 
00895                 // We have to initialize the guide particles during the first step.
00896                 cSMDInitGuides = 0;
00897                 
00898                 cSMDnumInstances = molPtr->get_cSMDnumInst();
00899                 cSMDindex = molPtr->get_cSMDindex();
00900                 cSMDindxLen = molPtr->get_cSMDindxLen();
00901                 cSMDpairs = molPtr->get_cSMDpairs(); 
00902                 cSMDKs = molPtr->get_cSMDKs();
00903                 cSMDVels = molPtr->get_cSMDVels();
00904                 cSMDcoffs = molPtr->get_cSMDcoffs();
00905                 
00906                 cSMDguides = new Position[cSMDnumInstances];
00907                 cSMDForces = new Force[cSMDnumInstances];
00908         }
00909                 
00910                 
00911         DebugM(4,"Initializing DCD file for charge information." << std::endl);
00912         
00913         // Initializes output DCD file for charge information.
00914         if (simParams->qmOutFreq > 0) {
00915             std::string dcdOutFileStr;
00916             dcdOutFileStr.clear();
00917             dcdOutFileStr.append(simParams->outputFilename) ;
00918             dcdOutFileStr.append(".qdcd") ;
00919             dcdOutFile = open_dcd_write(dcdOutFileStr.c_str()) ;
00920             
00921             if (dcdOutFile == DCD_FILEEXISTS) {
00922                 iout << iERROR << "DCD file " << dcdOutFile << " already exists!!\n" << endi;
00923                 NAMD_err("Could not write QM charge DCD file.");
00924             } else if (dcdOutFile < 0) {
00925                 iout << iERROR << "Couldn't open DCD file " << dcdOutFile << ".\n" << endi;
00926                 NAMD_err("Could not write QM charge DCD file.");
00927             } else if (! dcdOutFile) {
00928                 NAMD_bug("ComputeQMMgr::recvPartQM open_dcd_write returned fileid of zero");
00929             }
00930             
00931             timeStep = simParams->firstTimestep;
00932             
00933             int NSAVC, NFILE, NPRIV, NSTEP;
00934             NSAVC = simParams->qmOutFreq;
00935             NPRIV = timeStep;
00936             NSTEP = NPRIV - NSAVC;
00937             NFILE = 0;
00938             
00939             //  Write out the header
00940             int ret_code = write_dcdheader(dcdOutFile, dcdOutFileStr.c_str(),
00941             numQMAtms, NFILE, NPRIV, NSAVC, NSTEP,
00942             simParams->dt/TIMEFACTOR, 0);
00943 
00944             if (ret_code<0) {
00945                 NAMD_err("Writing of DCD header failed!!");
00946             }
00947             
00948             // The DCD write function takes 3 independent arrays for X, Y and Z
00949             // coordinates, but we allocate one and send in the pieces.
00950             outputData = new Real[3*numQMAtms];
00951         }
00952         
00953         DebugM(4,"Initializing DCD file for position information." << std::endl);
00954         // Initializes output DCD file for position information.
00955         if (simParams->qmPosOutFreq > 0) {
00956             std::string dcdPosOutFileStr;
00957             dcdPosOutFileStr.clear();
00958             dcdPosOutFileStr.append(simParams->outputFilename) ;
00959             dcdPosOutFileStr.append(".QMonly.dcd") ;
00960             dcdPosOutFile = open_dcd_write(dcdPosOutFileStr.c_str()) ;
00961             
00962             if (dcdPosOutFile == DCD_FILEEXISTS) {
00963                 iout << iERROR << "DCD file " << dcdPosOutFile << " already exists!!\n" << endi;
00964                 NAMD_err("Could not write QM charge DCD file.");
00965             } else if (dcdPosOutFile < 0) {
00966                 iout << iERROR << "Couldn't open DCD file " << dcdPosOutFile << ".\n" << endi;
00967                 NAMD_err("Could not write QM charge DCD file.");
00968             } else if (! dcdPosOutFile) {
00969                 NAMD_bug("ComputeQMMgr::recvPartQM open_dcd_write returned fileid of zero");
00970             }
00971             
00972             timeStep = simParams->firstTimestep;
00973             
00974             int NSAVC, NFILE, NPRIV, NSTEP;
00975             NSAVC = simParams->qmOutFreq;
00976             NPRIV = timeStep;
00977             NSTEP = NPRIV - NSAVC;
00978             NFILE = 0;
00979             
00980             //  Write out the header
00981             int ret_code = write_dcdheader(dcdPosOutFile, dcdPosOutFileStr.c_str(),
00982             numQMAtms, NFILE, NPRIV, NSAVC, NSTEP,
00983             simParams->dt/TIMEFACTOR, 0);
00984 
00985             if (ret_code<0) {
00986                 NAMD_err("Writing of DCD header failed!!");
00987             }
00988             
00989             // The DCD write function takes 3 independent arrays for X, Y and Z
00990             // coordinates, but we allocate one and send in the pieces.
00991             outputData = new Real[3*numQMAtms];
00992         }
00993         
00994         // Prepares list of PEs which will run the QM software
00995         int simsPerNode = simParams->qmSimsPerNode ;
00996         int zeroNodeSize = CmiNumPesOnPhysicalNode(0);
00997         
00998         // Check if the node has enought PEs to run the requested number of simulations.
00999         if ( zeroNodeSize < simsPerNode ) {
01000             iout << iERROR << "There are " << zeroNodeSize << " cores pernode, but "
01001             << simsPerNode << " QM simulations per node were requested.\n" << endi ;
01002             NAMD_die("Error preparing QM simulations.");
01003         }
01004         
01005         DebugM(4,"Preparing PE list for QM execution.\n");
01006         qmPEs.clear(); // Making sure its empty.
01007         
01008         int numNodes = CmiNumPhysicalNodes();
01009         int numPlacedQMGrps = 0;
01010         int placedOnNode = 0;
01011         int nodeIt = 0 ;
01012         
01013         // The default is to only run on rank zero.
01014         if ( simsPerNode <= 0 ) {
01015             qmPEs.push_back(0);
01016             numPlacedQMGrps = 1;
01017         }
01018         
01019         while ( (numPlacedQMGrps < numQMGrps) && (simsPerNode > 0) ) {
01020             
01021             // If we searched all nodes, break the loop.
01022             if (nodeIt == numNodes) {
01023                 break;
01024             }
01025             
01026             int *pelist;
01027             int nodeSize;
01028             CmiGetPesOnPhysicalNode(CmiPhysicalNodeID(nodeIt), &pelist, &nodeSize);
01029             
01030             DebugM(4,"Checking node " << nodeIt +1 << " out of " << numNodes
01031             << " (" << nodeSize << " PEs: " << pelist[0] << " to " 
01032             << pelist[nodeSize-1] << ")." << std::endl );
01033             
01034             for ( int i=0; i<nodeSize; ++i ) {
01035                 
01036                 // Check if the PE has patches. We only run on PEs with patches!
01037                 if ( (PatchMap::Object())->numPatchesOnNode(pelist[i]) == 0 ) {
01038                     DebugM(1,"PE " << pelist[i] << " has no patches! We'll skip it." 
01039                     << std::endl);
01040                     continue ;
01041                 }
01042                 
01043                 // Add the target PE on the target node to the list
01044                 // of PEs which will carry QM simulations.
01045                 qmPEs.push_back(pelist[i]);
01046                 
01047                 DebugM(1,"Added PE to QM execution list: " << pelist[i]  << "\n");
01048                 
01049                 numPlacedQMGrps++;
01050                 placedOnNode++;
01051                 
01052                 if (placedOnNode == simsPerNode) {
01053                     DebugM(1,"Placed enought simulations on this node.\n");
01054                     break;
01055                 }
01056                 
01057                 
01058             }
01059             
01060             nodeIt++ ;
01061             placedOnNode = 0;
01062         }
01063         
01064         if ( numPlacedQMGrps < numQMGrps ) {
01065             iout << iWARN << "Could not compute all QM groups in parallel.\n" << endi ;
01066         }
01067         
01068         iout << iINFO << "List of ranks running QM simulations: " << qmPEs[0] ;
01069         for (int i=1; i < qmPEs.size(); i++) {
01070             iout << ", " << qmPEs[i] ;
01071         }
01072         iout << ".\n" << endi;
01073         
01074     }
01075     
01076     DebugM(1,"Receiving from rank " << msg->sourceNode
01077     << " a total of " << msg->numAtoms << " QM atoms." << std::endl);
01078     
01079     // In case we are NOT using point charges but there are QM-MM bonds,
01080     // test each QM message for MM1 atoms.
01081     if (meNumMMIndx > 0) {
01082         
01083         ResizeArray< ComputeQMAtom > tmpAtm;
01084         ComputeQMAtom *data_ptr = msg->coord;
01085         
01086         for (int i=0; i<msg->numAtoms; i++) {
01087             if (data_ptr[i].vdwType < 0) {
01088                 tmpAtm.add(data_ptr[i]) ;
01089             }
01090         }
01091         
01092         QMPntChrgMsg *pntChrgMsg = new (tmpAtm.size(), 0) QMPntChrgMsg;
01093         pntChrgMsg->sourceNode = msg->sourceNode ;
01094         pntChrgMsg->numAtoms = tmpAtm.size() ;
01095         ComputeQMPntChrg* newPCData = pntChrgMsg->coord ;
01096         
01097         QMCoordMsg *newMsg = msg;
01098         
01099         if (tmpAtm.size() > 0) {
01100             
01101             newMsg = new (msg->numAtoms - tmpAtm.size(),0, 0) QMCoordMsg;
01102             newMsg->sourceNode = msg->sourceNode ;
01103             newMsg->numAtoms = msg->numAtoms - tmpAtm.size() ;
01104             newMsg->timestep = msg->timestep ;
01105             ComputeQMAtom *newMsgData = newMsg->coord;
01106             
01107             for (int i=0; i<msg->numAtoms; i++) {
01108                 if (data_ptr[i].vdwType < 0) {
01109                     newPCData->position = data_ptr[i].position ;
01110                     newPCData->charge = data_ptr[i].charge ;
01111                     newPCData->id = data_ptr[i].id ;
01112                     newPCData->qmGrpID = data_ptr[i].qmGrpID ;
01113                     newPCData->homeIndx = data_ptr[i].homeIndx ;
01114                     newPCData->dist = 0 ;
01115                     newPCData->mass = 0 ;
01116                     newPCData->vdwType = 0 ;
01117                     newPCData++;
01118                 }
01119                 else {
01120                     *newMsgData = data_ptr[i] ;
01121                     newMsgData++;
01122                 }
01123             }
01124             
01125             delete msg;
01126             
01127         }
01128         
01129         qmCoordMsgs[numArrivedQMMsg] = newMsg;
01130         ++numArrivedQMMsg;
01131         
01132         pntChrgCoordMsgs[numArrivedPntChrgMsg] = pntChrgMsg;
01133         ++numArrivedPntChrgMsg;
01134     }
01135     else {
01136         qmCoordMsgs[numArrivedQMMsg] = msg;
01137         ++numArrivedQMMsg;
01138     }
01139     
01140     if ( numArrivedQMMsg < numSources ) 
01141         return;
01142     
01143     // Now that all messages arrived, get all QM positions.
01144     for (int msgIt=0; msgIt < numArrivedQMMsg; msgIt++){
01145         
01146         DebugM(1, "Getting positions for " << qmCoordMsgs[msgIt]->numAtoms 
01147         << " QM atoms in this message." << std::endl);
01148         
01149         for ( int i=0; i < qmCoordMsgs[msgIt]->numAtoms; ++i ) {
01150             qmCoord[qmAtmIter] = qmCoordMsgs[msgIt]->coord[i];
01151             qmAtmIter++;
01152         }
01153     }
01154     
01155     if (qmAtmIter != numQMAtms) {
01156         iout << iERROR << "The number of QM atoms received (" << qmAtmIter
01157         << ") is different than expected: " << numQMAtms << "\n" << endi;
01158         NAMD_err("Problems broadcasting QM atoms.");
01159     }
01160     
01161     // Resets the counter for the next step.
01162     numArrivedQMMsg = 0;
01163     qmAtmIter = 0;
01164     
01165     timeStep = qmCoordMsgs[0]->timestep;
01166     
01167     // Makes sure there is no trash or old info in the force array.
01168     // This is where we accumulate forces from the QM software and our own
01169     // Coulomb forces. It will have info on QM atoms and point charges only.
01170     for (int i=0; i<numAtoms; ++i ) {
01171         force[i].force = 0;  // this assigns 0 (ZERO) to all componenets of the vector.
01172         force[i].replace = 0;  // By default, no atom has all forces replaced.
01173         force[i].homeIndx = -1;  // To prevent errors from sliping.
01174         force[i].charge = 0;
01175         force[i].id = i;        // Initializes the variable since not all forces are sent to all ranks.
01176     }
01177     
01178     
01179     for (int i=0; i<numQMAtms; i++) {
01180         // Each force receives the home index of its atom with respect to the 
01181         // local set of atoms in each node.
01182         if (force[qmCoord[i].id].homeIndx != -1
01183             && force[qmCoord[i].id].homeIndx != qmCoord[i].homeIndx
01184         ) {
01185             iout << iERROR << "Overloading QM atom " 
01186             << qmCoord[i].id << "; home index: " 
01187             << force[qmCoord[i].id].homeIndx << " with " << qmCoord[i].homeIndx
01188             << "\n" << endi ;
01189             NAMD_die("Error preparing QM simulations.");
01190         }
01191         
01192         force[qmCoord[i].id].homeIndx = qmCoord[i].homeIndx;
01193         // Each force on QM atoms has an indicator so NAMD knows if all 
01194         // NAMD forces should be erased and completely overwritten by the
01195         // external QM forces.
01196         force[qmCoord[i].id].replace = replaceForces;
01197     }
01198     
01199     if (noPC) {
01200         // this pointer should be freed in rank zero, after receiving it.
01201         QMPntChrgMsg *pntChrgMsg = new (0, 0) QMPntChrgMsg;
01202         pntChrgMsg->sourceNode = CkMyPe();
01203         pntChrgMsg->numAtoms = 0;
01204         
01205         QMProxy[0].recvPntChrg(pntChrgMsg);
01206     }
01207     else {
01208         // The default mode is to update the poitn charge selection at every step.
01209         int pcSelMode = PCMODEUPDATESEL;
01210         
01211         int msgPCListSize = 0;
01212         // We check wether point charges are to be selected with a stride.
01213         if ( qmPCFreq > 0 ) {
01214             if (timeStep % qmPCFreq == 0) {
01215                 // Marks that the PC list determined in this step will
01216                 // be broadcast on the *next* step, and kept for the following
01217                 // qmPCFreq-1 steps.
01218                 resendPCList = true;
01219                 
01220                 // Clears the list since we don't know how many charges 
01221                 // will be selected this time.
01222                 delete [] pcIDList;
01223             }
01224             else {
01225                 // If the PC selection is not to be updated in this step,
01226                 // inform the nodes that the previous list (or the updated
01227                 // list being sent in this message) is to be used and only 
01228                 // updated positions will be returned.
01229                 pcSelMode = PCMODEUPDATEPOS;
01230                 
01231                 // If this is the first step after a PC re-selection, all 
01232                 // ranks receive the total list, since atoms may move between
01233                 // PEs in between PC re-assignments (if they are far enought apart
01234                 // or if you are unlucky)
01235                 if (resendPCList) {
01236                     msgPCListSize = pcIDListSize;
01237                     resendPCList = false;
01238                 }
01239             }
01240         }
01241         
01242         // In case we are using custom selection of point charges, indicate the mode.
01243         if (simParams->qmCustomPCSel)
01244             pcSelMode = PCMODECUSTOMSEL;
01245         
01246         DebugM(1,"Broadcasting current positions of QM atoms.\n");
01247         for ( int j=0; j < numSources; ++j ) {
01248             // This pointer will be freed in the receiving rank.
01249             QMCoordMsg *qmFullMsg = new (numQMAtms, msgPCListSize, 0) QMCoordMsg;
01250             qmFullMsg->sourceNode = CkMyPe();
01251             qmFullMsg->numAtoms = numQMAtms;
01252             qmFullMsg->pcSelMode = pcSelMode;
01253             qmFullMsg->numPCIndxs =  msgPCListSize;
01254             ComputeQMAtom *data_ptr = qmFullMsg->coord;
01255             int *msgPCListP = qmFullMsg->pcIndxs;
01256             
01257             for (int i=0; i<numQMAtms; i++) {
01258                 data_ptr->position = qmCoord[i].position;
01259                 data_ptr->charge = qmCoord[i].charge;
01260                 data_ptr->id = qmCoord[i].id;
01261                 data_ptr->qmGrpID = qmCoord[i].qmGrpID;
01262                 data_ptr->homeIndx = -1; // So there is no mistake that there is no info here.
01263                 data_ptr++;
01264             }
01265             
01266             for (int i=0; i<msgPCListSize; i++) {
01267                 msgPCListP[i] = pcIDList[i] ;
01268             }
01269             
01270             #ifdef DEBUG_QM
01271             if (j == 0)
01272                 Write_PDB("qmMsg.pdb", qmFullMsg) ;
01273             #endif
01274             
01275             // The messages are deleted later, we will need them.
01276             QMProxy[qmCoordMsgs[j]->sourceNode].recvFullQM(qmFullMsg);
01277         }
01278     }
01279 }

void ComputeQMMgr::recvPntChrg ( QMPntChrgMsg msg  ) 

Definition at line 1910 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, QMAtomData::charge, QMGrpCalcMsg::charge, ComputeQMAtom::charge, ResizeArray< Elem >::clear(), QMGrpCalcMsg::configLines, Node::configList, QMGrpCalcMsg::constants, COULOMB, custom_ComputeQMAtom_Less(), SimParameters::cutoff, QMGrpCalcMsg::cutoff, QMGrpCalcMsg::data, StringList::data, DebugM, SimParameters::dielectric, QMAtomData::dist, QMAtomData::element, endi(), QMGrpCalcMsg::execPath, ConfigList::find(), SortedArray< Elem >::find(), SimParameters::firstTimestep, Molecule::get_qmDummyBondVal(), Molecule::get_qmDummyElement(), Molecule::get_qmElements(), Molecule::get_qmGrpBonds(), Molecule::get_qmGrpNumBonds(), Molecule::get_qmGrpSizes(), Molecule::get_qmMMBond(), Molecule::get_qmMMBondedIndx(), Molecule::get_qmMMChargeTarget(), Molecule::get_qmMMNumTargs(), QMGrpCalcMsg::grpIndx, QMForce::homeIndx, QMAtomData::id, iERROR(), iout, SortedArray< Elem >::load(), QMGrpCalcMsg::multiplicity, NAMD_die(), StringList::next, SimParameters::nonbondedScaling, QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMPntChrgMsg::numAtoms, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, Node::Object(), QMGrpCalcMsg::pcScheme, QMGrpCalcMsg::peIter, SimParameters::PMEEwaldCoefficient, QMGrpCalcMsg::PMEEwaldCoefficient, SimParameters::PMEOn, QMGrpCalcMsg::PMEOn, QMAtomData::position, QMGrpCalcMsg::prepProc, QMGrpCalcMsg::prepProcOn, QMGrpCalcMsg::qmAtmChrgMode, QMATOMTYPE_DUMMY, QMATOMTYPE_QM, SimParameters::qmBaseDir, SimParameters::qmChrgMode, SimParameters::qmExecPath, SimParameters::qmFormat, QMFormatMOPAC, QMFormatORCA, QMFormatUSR, QMLENTYPE, SimParameters::qmMOPACAddConfigChrg, SimParameters::qmPCScheme, SimParameters::qmPCSwitchOn, SimParameters::qmPCSwitchType, QMPCTYPE_CLASSICAL, QMPCTYPE_EXTRA, QMPCTYPE_IGNORE, SimParameters::qmPrepProc, SimParameters::qmPrepProcOn, QMRATIOTYPE, SimParameters::qmSecProc, SimParameters::qmSecProcOn, QMGrpCalcMsg::secProc, QMGrpCalcMsg::secProcOn, ResizeArray< Elem >::size(), sort, SortedArray< Elem >::sort(), QMGrpCalcMsg::swdist, QMGrpCalcMsg::switching, SimParameters::switchingDist, QMGrpCalcMsg::switchType, QMGrpCalcMsg::timestep, QMAtomData::type, and Vector::unit().

01910                                                 {
01911     
01912     // All the preparation that used to be in this function was moved to 
01913     // recvPartQM, which is called first in rank zero.
01914     
01915     if (noPC) {
01916         // Even if we have QM-MM bonds, the point charge messages 
01917         // are handled in recvPartQM
01918         delete msg;
01919     }
01920     else {
01921         pntChrgCoordMsgs[numArrivedPntChrgMsg] = msg;
01922         ++numArrivedPntChrgMsg;
01923 
01924         if ( numArrivedPntChrgMsg < numSources ) 
01925             return;
01926     }
01927     
01928     // Resets counters for next step.
01929     numRecQMRes = 0;
01930     
01931     totalEnergy = 0;
01932     
01933     for ( int k=0; k<3; ++k )
01934         for ( int l=0; l<3; ++l )
01935             totVirial[k][l] = 0;
01936     
01937     // ALL DATA ARRIVED --- CALCULATE FORCES
01938     
01939     const char *const *const elementArray = molPtr->get_qmElements() ;
01940     const char *const *const dummyElmArray = molPtr->get_qmDummyElement();
01941     const int *const qmGrpSize = molPtr->get_qmGrpSizes();
01942     
01943     const BigReal *const dummyBondVal = molPtr->get_qmDummyBondVal();
01944     const int *const grpNumBonds = molPtr->get_qmGrpNumBonds() ;
01945     const int *const *const qmMMBond = molPtr->get_qmMMBond() ;
01946     const int *const *const qmGrpBonds = molPtr->get_qmGrpBonds() ;
01947     const int *const *const qmMMBondedIndx = molPtr->get_qmMMBondedIndx() ;
01948     const int *const *const chargeTarget = molPtr->get_qmMMChargeTarget() ;
01949     const int *const numTargs = molPtr->get_qmMMNumTargs() ;
01950     
01951     BigReal constants = COULOMB * simParams->nonbondedScaling / (simParams->dielectric) ;
01952     // COULOMB is in kcal*Angs/(mol*e^2)
01953 //     BigReal constants = COULOMB ;
01954     
01955     if ( qmPCFreq > 0 ) {
01956         DebugM(4,"Using point charge stride of " << qmPCFreq << "\n")
01957         // In case we are skiping steps between PC selection, store a main list
01958         // in rank zero for future steps. Then we only update positions untill
01959         // we reach a new "qmPCFreq" step, when a new PC selection is made.
01960         
01961         if (timeStep % qmPCFreq == 0) {
01962             DebugM(4,"Loading a new selection of PCs.\n")
01963             
01964             // We only re-set this arrya in a step where a new PC selection is made.
01965             pntChrgMsgVec.clear();
01966             for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
01967                 // Accumulates the message point charges into a local vector.
01968                 for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
01969                     pntChrgMsgVec.push_back(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
01970                 }
01971             }
01972             
01973             // This fast array is created to send the entire PC IDs list to all
01974             // PEs in the next step.
01975             pcIDListSize = pntChrgMsgVec.size();
01976             pcIDList = new int[pcIDListSize] ;
01977             for (size_t i=0; i<pcIDListSize; i++) {
01978                 pcIDList[i] = pntChrgMsgVec[i].id;
01979                 
01980                 // Loads home indexes of MM atoms received as point charges.
01981                 // Each force receives the home index of its atom with respect to the 
01982                 // local set of atoms in each node.
01983                 force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
01984             }
01985         }
01986         else {
01987             DebugM(4,"Updating position/homeIdex of old PC selection.\n")
01988             
01989             // We build a sorted array so that we can quickly access the unordered
01990             // data we just received, and update positions of the same selection
01991             // of point charges.
01992             pcDataSort.clear();
01993             for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
01994                 // Accumulates the message point charges into a local sorted array.
01995                 for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
01996                     pcDataSort.load(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
01997                 }
01998             }
01999             pcDataSort.sort();
02000             
02001             iout << "Loaded new positions in sorted array: " << pcDataSort.size() << "\n" << endi;
02002             
02003             // If we are using a set of point charges that was selected in a
02004             // previous step, we update the PC POSITION ONLY.
02005             for (size_t i=0; i<pntChrgMsgVec.size(); i++) {
02006                 
02007                 auto pntr = pcDataSort.find(pntChrgMsgVec[i]);
02008                 
02009                 pntChrgMsgVec[i].position = pntr->position ;
02010                 pntChrgMsgVec[i].homeIndx = pntr->homeIndx ;
02011                 
02012                 // Loads home indexes of MM atoms received as point charges.
02013                 // Each force receives the home index of its atom with respect to the 
02014                 // local set of atoms in each node.
02015                 force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
02016             }
02017         }
02018     }
02019     else {
02020         DebugM(4,"Updating PCs at every step.\n")
02021         
02022         pntChrgMsgVec.clear();
02023         for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
02024             // Accumulates the message point charges into a local vector.
02025             for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
02026                 pntChrgMsgVec.push_back(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
02027             }
02028         }
02029         
02030         // Loads home indexes of MM atoms received as point charges.
02031         for (size_t i=0; i<pntChrgMsgVec.size(); i++) {
02032             // Each force receives the home index of its atom with respect to the 
02033             // local set of atoms in each node.
02034             
02035             #ifdef DEBUG_QM
02036             if (force[pntChrgMsgVec[i].id].homeIndx != -1 
02037                 and force[pntChrgMsgVec[i].id].homeIndx != pntChrgMsgVec[i].homeIndx
02038             ) {
02039                 iout << iERROR << "Overloading point charge " 
02040                 << pntChrgMsgVec[i].id << "; home index: " 
02041                 << force[pntChrgMsgVec[i].id].homeIndx << " with " << pntChrgMsgVec[i].homeIndx
02042                 << "\n" << endi ;
02043                 NAMD_die("Error preparing QM simulations.");
02044             }
02045             #endif
02046             
02047             force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
02048         }
02049     }
02050     
02051     // Reset counter for next timestep
02052     numArrivedPntChrgMsg = 0;
02053     
02054     DebugM(4,"A total of " << pntChrgMsgVec.size() 
02055     << " point charges arrived." << std::endl);
02056     
02057     DebugM(4,"Starting QM groups processing.\n");
02058     
02059     QMAtmVec grpQMAtmVec;
02060     QMPCVec grpPntChrgVec;
02061     
02062     // Final set of charges, created or not, that is sent to the QM software.
02063     // This set will depend on how QM-MM bonds are processed and presented to the
02064     // QM region.
02065     QMPCVec grpAppldChrgVec;
02066     
02067     // Vector of dummy atoms created to treat QM-MM bonds.
02068     std::vector<dummyData> dummyAtoms ;
02069 
02070     // Initializes the loop for receiving the QM results.
02071     thisProxy[0].recvQMResLoop() ;
02072     
02073     // Iterator for target PE where QM simulations will run.
02074     int peIter = 0;
02075     
02076     for (int grpIter = 0; grpIter < numQMGrps; grpIter++) {
02077         
02078         grpQMAtmVec.clear();
02079         grpPntChrgVec.clear();
02080         grpAppldChrgVec.clear();
02081         dummyAtoms.clear();
02082         
02083         DebugM(4,"Calculating QM group " << grpIter +1 
02084         << " (ID: " << grpID[grpIter] << ")." << std::endl);
02085         
02086         DebugM(4,"Compiling Config Lines into one string for message...\n");
02087         
02088         // This will hold a big sting with all configuration lines the user supplied.
02089         int lineIter = 0 ;
02090         std::string configLines ;
02091         StringList *current = Node::Object()->configList->find("QMConfigLine");
02092         for ( ; current; current = current->next ) {
02093             std::string confLineStr(current->data);
02094             
02095             // In case we need to add charges to MOPAC command line.
02096             if (simParams->qmFormat == QMFormatMOPAC && simParams->qmMOPACAddConfigChrg && lineIter == 0) {
02097                 std::ostringstream itosConv ;
02098                 itosConv << grpChrg[grpIter] ;
02099                 confLineStr.append( " CHARGE=" );
02100                 confLineStr.append( itosConv.str() );
02101                 
02102             }
02103             
02104             configLines.append(confLineStr);
02105             configLines.append("\n");
02106             
02107             lineIter++;
02108         }
02109         
02110         DebugM(4,"Determining point charges...\n");
02111         
02112         Real qmTotalCharge = 0;
02113         // Loads the QM atoms for this group.
02114         for (int qmIt=0; qmIt<numQMAtms; qmIt++){
02115             if (qmCoord[qmIt].qmGrpID == grpID[grpIter]) {
02116                 grpQMAtmVec.push_back(qmCoord[qmIt]);
02117                 qmTotalCharge += qmCoord[qmIt].charge;
02118             }
02119         }
02120         if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
02121             qmTotalCharge = roundf(qmTotalCharge) ;
02122         }
02123         
02124         // Sorts the vector so that the QM message is always built with the same order of atoms.
02125         std::sort(grpQMAtmVec.begin(), grpQMAtmVec.end(), custom_ComputeQMAtom_Less);
02126         
02127         Real pcTotalCharge = 0;
02128         // Loads the point charges to a local vector for this QM group.
02129         for (auto pntChrgIt = pntChrgMsgVec.begin(); 
02130              pntChrgIt != pntChrgMsgVec.end(); pntChrgIt++) {
02131             if ((*pntChrgIt).qmGrpID == grpID[grpIter] ) {
02132                 grpPntChrgVec.push_back( (*pntChrgIt) );
02133                 pcTotalCharge += (*pntChrgIt).charge;
02134             }
02135         }
02136         if ((fabsf(roundf(pcTotalCharge) - pcTotalCharge) <= 0.001f) ) {
02137             pcTotalCharge = roundf(pcTotalCharge) ;
02138         }
02139         
02140         #ifdef DEBUG_QM
02141         if (grpQMAtmVec.size() != qmGrpSize[grpIter]) {
02142             iout << iERROR << "The number of QM atoms received for group " << grpID[grpIter]
02143             << " does not match the expected: " << grpQMAtmVec.size()
02144             << " vs " << qmGrpSize[grpIter] << "\n" << endi ;
02145             
02146             NAMD_die("Error processing QM group.");
02147         }
02148         #endif
02149         
02150         DebugM(1,"Found " << grpPntChrgVec.size() << " point charges for QM group " 
02151         << grpIter << " (ID: " << grpID[grpIter] << "; Num QM atoms: " 
02152         << grpQMAtmVec.size() <<  "; Num QM-MM bonds: " 
02153         << grpNumBonds[grpIter] << ")" << std::endl);
02154         
02155         DebugM(1,"Total QM charge: " << qmTotalCharge 
02156         << "; Total point-charge charge: " << pcTotalCharge << std::endl);
02157         
02158         // If we have a frequency for LSS update, check if we shoudl do it in 
02159         // the current time step.
02160         if ( qmLSSFreq > 0 && ((timeStep + 1) % qmLSSFreq == 0 )) {
02161             lssUpdate(grpIter, grpQMAtmVec, grpPntChrgVec);
02162         }
02163         
02164         // This function checks data and treats the charge (and existence) of
02165         // the MM atoms in and around QM-MM bonds. It is only executed in 
02166         // electrostatic embeding QM/MM simulations.
02167         if (! noPC ) {
02168             procBonds(grpNumBonds[grpIter], qmGrpBonds[grpIter], 
02169                      qmMMBondedIndx[grpIter], 
02170                      chargeTarget, numTargs, 
02171                      grpPntChrgVec, grpAppldChrgVec) ;
02172         }
02173         else {
02174             grpAppldChrgVec = grpPntChrgVec;
02175         }
02176         
02177         // For future use, we get the pairs of indexes of QM atoms and MM atoms which are
02178         // bound in QM-MM bonds.
02179         std::vector< std::pair<int,int> > qmPCLocalIndxPairs ;
02180         
02181         // Create and position dummy atoms.
02182         Position mmPos(0), qmPos(0);
02183         for (int dummyIt = 0; dummyIt < grpNumBonds[grpIter]; dummyIt++) {
02184             
02185             int qmMMBondIndx = qmGrpBonds[grpIter][dummyIt] ;
02186             
02187             BigReal bondVal = dummyBondVal[qmMMBondIndx] ;
02188             
02189             int mmAtmIndx = qmMMBond[qmMMBondIndx][0] ;
02190             int qmAtmIndx = qmMMBond[qmMMBondIndx][1] ;
02191             
02192             // Sicne we don't know in which patch/node the QM atom is, or the 
02193             // order in which they will arrive in rank zero, we have
02194             // no direct index to it.
02195             #ifdef DEBUG_QM
02196             bool missingQM = true, missingMM = true;
02197             #endif
02198             size_t qmIt ;
02199             for (qmIt=0; qmIt<grpQMAtmVec.size(); qmIt++){
02200                 if (grpQMAtmVec[qmIt].id == qmAtmIndx) {
02201                     qmPos = grpQMAtmVec[qmIt].position;
02202                     
02203                     #ifdef DEBUG_QM
02204                     missingQM = false;
02205                     #endif
02206                     
02207                     break;
02208                 }
02209             }
02210             // The same is true about the MM atom to which the QM atom is bound,
02211             // we must look
02212             size_t pcIt;
02213             for (pcIt=0; pcIt < grpPntChrgVec.size(); pcIt++) {
02214                 if (grpPntChrgVec[pcIt].id == mmAtmIndx) {
02215                     mmPos = grpPntChrgVec[pcIt].position ;
02216                     
02217                     #ifdef DEBUG_QM
02218                     missingMM = false;
02219                     #endif
02220                     
02221                     break;
02222                 }
02223             }
02224             
02225             qmPCLocalIndxPairs.push_back(std::pair<int,int>(qmIt, pcIt) );
02226             
02227             #ifdef DEBUG_QM
02228             // Checks if the MM atom was included as a 
02229             // point charge due proximity to a QM atom, and if the QM atom arrived.
02230             if ( missingMM or missingQM ) {
02231                 // If it wasn't, there is an error somwhere.
02232                 
02233                 if (missingMM) {
02234                     iout << iERROR << "The MM atom " << mmAtmIndx
02235                     << " is bound to a QM atom, but it was not selected as a poitn charge."
02236                     << " Check your cutoff radius!\n" << endi ;
02237                     
02238                     NAMD_die("Error in QM-MM bond processing.");
02239                 }
02240                 if (missingQM) {
02241                     iout << iERROR << "The QM atom " << qmAtmIndx
02242                     << " is bound to an MM atom, but it was not sent to rank zero for processing."
02243                     << " Check your configuration!\n" << endi ;
02244                     
02245                     NAMD_die("Error in QM-MM bond processing.");
02246                 }
02247             }
02248             #endif
02249             
02250             Vector bondVec = mmPos - qmPos ;
02251             
02252             if (bondValType == QMLENTYPE) {
02253                 // If a length is defined by the user, or a default len
02254                 // is used, we calculate the unit vector for the displacement
02255                 // and multiply by the desired length in order 
02256                 // to get the final dummy atom position relative to the
02257                 // QM atom.
02258                 bondVec = bondVec.unit() ;
02259                 bondVec *= bondVal ;
02260             }
02261             else if (bondValType == QMRATIOTYPE) {
02262                 // If distance a ratio was defined by the user, then
02263                 // the displacement vector is multiplied by that ratio
02264                 // to get the final dummy atom position relative to the
02265                 // QM atom.
02266                 bondVec *= bondVal ;
02267             }
02268             
02269             Position dummyPos = qmPos + bondVec;
02270             
02271             DebugM(1,"Creating dummy atom " << dummyPos << " ; QM ind: " 
02272             << qmIt << " ; PC ind: " << pcIt << std::endl);
02273             
02274             dummyAtoms.push_back(dummyData(dummyPos, qmIt, qmMMBondIndx)) ;
02275             
02276         }
02277         
02278         DebugM(3, "Creating data for " << grpQMAtmVec.size() << " QM atoms " 
02279         << dummyAtoms.size() << " dummy atoms " << grpPntChrgVec.size()
02280         << " real point charges and " << grpAppldChrgVec.size() - grpPntChrgVec.size()
02281         << " virtual point charges\n") ;
02282         
02283         int dataSize = grpQMAtmVec.size() + dummyAtoms.size() + grpAppldChrgVec.size();
02284         QMGrpCalcMsg *msg = new (dataSize, configLines.size(), 0) QMGrpCalcMsg;
02285         msg->timestep = timeStep;
02286         msg->grpIndx = grpIter;
02287         msg->peIter = peIter;
02288         msg->charge = grpChrg[grpIter];
02289         msg->multiplicity = grpMult[grpIter];
02290         msg->numQMAtoms = grpQMAtmVec.size();
02291         msg->numAllAtoms = grpQMAtmVec.size() + dummyAtoms.size();
02292         msg->numRealPntChrgs = grpPntChrgVec.size(); // The original set of classical atoms.
02293         msg->numAllPntChrgs = grpAppldChrgVec.size(); // The extended set with virtual point charges.
02294         msg->secProcOn = simParams->qmSecProcOn ;
02295         msg->constants = constants;
02296         msg->PMEOn = simParams->PMEOn ;
02297         if (msg->PMEOn)
02298             msg->PMEEwaldCoefficient = simParams->PMEEwaldCoefficient ;
02299         msg->switching = simParams->qmPCSwitchOn;
02300         msg->switchType = simParams->qmPCSwitchType;
02301         msg->cutoff = simParams->cutoff;
02302         msg->swdist = simParams->switchingDist;
02303         msg->pcScheme = simParams->qmPCScheme;
02304         msg->qmAtmChrgMode = simParams->qmChrgMode;
02305         
02306         strncpy(msg->baseDir, simParams->qmBaseDir, 256);
02307         strncpy(msg->execPath, simParams->qmExecPath, 256);
02308         if (msg->secProcOn)
02309             strncpy(msg->secProc, simParams->qmSecProc, 256);
02310         
02311         if (simParams->qmPrepProcOn && (timeStep == simParams->firstTimestep)) {
02312             msg->prepProcOn = true;
02313             strncpy(msg->prepProc, simParams->qmPrepProc, 256);
02314         } else
02315             msg->prepProcOn = false;
02316         
02317         QMAtomData *dataP = msg->data;
02318         
02319         for (int i=0; i<grpQMAtmVec.size(); i++) {
02320             dataP->position = grpQMAtmVec[i].position ;
02321             dataP->charge = grpQMAtmVec[i].charge ;
02322             dataP->id = grpQMAtmVec[i].id ;
02323             dataP->bountToIndx = -1;
02324             dataP->type = QMATOMTYPE_QM ;
02325             strncpy(dataP->element,elementArray[grpQMAtmVec[i].id],3);
02326             dataP++;
02327         }
02328         
02329         for (int i=0; i<dummyAtoms.size(); i++) {
02330             dataP->position = dummyAtoms[i].pos ;
02331             dataP->charge = 0 ;
02332             dataP->id = -1 ;
02333             dataP->bountToIndx = dummyAtoms[i].qmInd ;
02334             dataP->type = QMATOMTYPE_DUMMY ;
02335             strncpy(dataP->element,dummyElmArray[dummyAtoms[i].bondIndx],3);
02336             dataP++;
02337         }
02338         
02339         for (int i=0; i<grpAppldChrgVec.size(); i++) {
02340             dataP->position = grpAppldChrgVec[i].position ;
02341             dataP->charge = grpAppldChrgVec[i].charge ;
02342             // Point charges created to handle QM-MM bonds will have an id of -1.
02343             dataP->id = grpAppldChrgVec[i].id ;
02344             dataP->bountToIndx = -1;
02345             dataP->dist = grpAppldChrgVec[i].dist ;
02346             // If we are loading the classical atoms' charges
02347             // the point charge type is 1, unless it is from an 
02348             // atom which is bound to a QM atom.
02349             if (i < grpPntChrgVec.size()) {
02350                 if (grpAppldChrgVec[i].qmGrpID == -1) {
02351                     dataP->type = QMPCTYPE_IGNORE ;
02352                 }
02353                 else if (grpAppldChrgVec[i].qmGrpID == -2) {
02354                     dataP->type = QMPCTYPE_CLASSICAL ;
02355                     dataP->bountToIndx = grpAppldChrgVec[i].homeIndx;
02356                 }
02357                 else {
02358                     dataP->type = QMPCTYPE_CLASSICAL ;
02359                 }
02360             }
02361             else {
02362                 // Extra charges are created to handle QM-MM bonds (if they exist).
02363                 dataP->type = QMPCTYPE_EXTRA ;
02364                 dataP->bountToIndx = grpAppldChrgVec[i].id;
02365             }
02366             dataP++;
02367         }
02368         
02369         QMAtomData *qmP = msg->data ;
02370         QMAtomData *pcP = msg->data + msg->numAllAtoms ;
02371         
02372         // With this, every QM atom knows to which MM atom is is bound, 
02373         // and vice versa. This will be usefull later on to prevent them from
02374         // feeling eachother's electrostatic charges AND to place the dummy
02375         // atom forces on the "real" atoms that form the bond.
02376         for( auto vecPtr  = qmPCLocalIndxPairs.begin(); 
02377                   vecPtr != qmPCLocalIndxPairs.end(); 
02378                   vecPtr++ ) {
02379             
02380             int qmLocInd = (*vecPtr).first;
02381             int pcLocInd = (*vecPtr).second;
02382             
02383             qmP[qmLocInd].bountToIndx = pcLocInd ;
02384             pcP[pcLocInd].bountToIndx = qmLocInd ;
02385         }
02386         
02387         
02388         strcpy(msg->configLines, configLines.c_str());
02389         
02390         if (cSMDon)
02391             calcCSMD(grpIter, msg->numQMAtoms, qmP, cSMDForces) ;
02392                 
02393         int targetPE = qmPEs[peIter] ;
02394         
02395         DebugM(4,"Sending QM group " << grpIter << " (ID " << grpID[grpIter] 
02396         << ") to PE " << targetPE << std::endl);
02397         
02398         switch (simParams->qmFormat) {
02399             // Creates the command line that will be executed according to the 
02400             // chosen QM software, as well as the input file with coordinates.
02401             case QMFormatORCA:
02402                 QMProxy[targetPE].calcORCA(msg) ;
02403                 break;
02404             
02405             case QMFormatMOPAC:
02406                 QMProxy[targetPE].calcMOPAC(msg) ;
02407                 break;
02408             
02409             case QMFormatUSR:
02410                 QMProxy[targetPE].calcUSR(msg) ;
02411                 break;
02412         }
02413         
02414         peIter++;
02415         
02416         if (peIter == qmPEs.size())
02417             peIter = 0;
02418     }
02419     
02420 }

void ComputeQMMgr::setCompute ( ComputeQM c  )  [inline]

Definition at line 406 of file ComputeQM.C.

00406 { QMCompute = c; }

void ComputeQMMgr::storeQMRes ( QMGrpResMsg resMsg  ) 

Definition at line 2422 of file ComputeQM.C.

References QMForce::charge, DebugM, endi(), QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, for(), QMForce::force, QMGrpResMsg::force, FORMAT(), QMGrpResMsg::grpIndx, iout, QMGrpResMsg::numForces, SimParameters::qmEnergyOutFreq, QMETITLE(), and QMGrpResMsg::virial.

02422                                                  {
02423     
02424 //     iout << iINFO << "Storing QM results for region " << resMsg->grpIndx  
02425 //     << " (ID: "  << grpID[resMsg->grpIndx] 
02426 //     << ") with original energy: " << endi;
02427 //     std::cout << std::fixed << std::setprecision(6) << resMsg->energyOrig << endi;
02428 //     iout << " Kcal/mol\n" << endi;
02429     
02430 //     if (resMsg->energyCorr != resMsg->energyOrig) {
02431 //         iout << iINFO << "PME corrected energy: " << endi;
02432 //         std::cout << std::fixed << std::setprecision(6) << resMsg->energyCorr << endi;
02433 //         iout << " Kcal/mol\n" << endi;
02434 //     }
02435     
02436     if ( (timeStep % simParams->qmEnergyOutFreq ) == 0 ) {
02437         iout << QMETITLE(timeStep);
02438         iout << FORMAT(grpID[resMsg->grpIndx]);
02439         iout << FORMAT(resMsg->energyOrig);
02440         if (resMsg->energyCorr != resMsg->energyOrig) iout << FORMAT(resMsg->energyCorr);
02441         iout << "\n\n" << endi;
02442     }
02443     
02444     totalEnergy += resMsg->energyCorr ;
02445     
02446     for ( int k=0; k<3; ++k )
02447         for ( int l=0; l<3; ++l )
02448             totVirial[k][l] += resMsg->virial[k][l];
02449     
02450     QMForce *fres = resMsg->force ;
02451     Real qmTotalCharge = 0;
02452     
02453     for (int i=0; i<resMsg->numForces; i++) {
02454         
02455         force[ fres[i].id ].force += fres[i].force;
02456         
02457         // Indicates the result is a QM atom, and we should get its updated charge.
02458         if (fres[i].replace == 1) {
02459             force[ fres[i].id ].charge =  fres[i].charge;
02460             qmTotalCharge += fres[i].charge;
02461         }
02462     }
02463     
02464     if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
02465         qmTotalCharge = roundf(qmTotalCharge) ;
02466     }
02467     
02468     // In case we are calculating cSMD forces, apply them now.
02469     if (cSMDon) {
02470         for ( int i=0; i < cSMDindxLen[resMsg->grpIndx]; i++ ) {
02471             int cSMDit = cSMDindex[resMsg->grpIndx][i];
02472             force[ cSMDpairs[cSMDit][0] ].force += cSMDForces[cSMDit];
02473         }
02474     }
02475 
02476     DebugM(4,"QM total charge received is " << qmTotalCharge << std::endl);
02477     
02478     DebugM(4,"Current accumulated energy is " << totalEnergy << std::endl);
02479     
02480     numRecQMRes++;
02481     
02482     delete resMsg;
02483 }


The documentation for this class was generated from the following file:

Generated on 21 Sep 2020 for NAMD by  doxygen 1.6.1