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  ) 

Definition at line 2783 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::charge, QMGrpCalcMsg::configLines, QMGrpCalcMsg::constants, QMGrpCalcMsg::data, DebugM, QMAtomData::element, endi(), QMGrpCalcMsg::execPath, iERROR(), if(), iout, j, NAMD_die(), NAMD_err(), QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, QMGrpCalcMsg::peIter, QMGrpCalcMsg::PMEOn, QMAtomData::position, QMPCTYPE_IGNORE, QMGrpCalcMsg::switching, QMAtomData::type, Vector::x, x, Vector::y, y, Vector::z, and z.

02784 {
02785     
02786     FILE *inputFile,*outputFile,*chrgFile;
02787     int iret;
02788     
02789     const size_t lineLen = 256;
02790     char *line = new char[lineLen];
02791     
02792     std::string qmCommand, inputFileName, outputFileName, pntChrgFileName;
02793     
02794     // For coulomb calculation
02795     BigReal constants = msg->constants;
02796     
02797     double gradient[3];
02798     
02799     DebugM(4,"Running MOPAC on PE " << CkMyPe() << std::endl);
02800     
02801     QMAtomData *atmP = msg->data ;
02802     QMAtomData *pcP = msg->data + msg->numAllAtoms ;
02803     
02804     // We create a pointer for PME correction, which may point to
02805     // a copy of the original point charge array, with unchanged charges,
02806     // or points to the original array in case no switching or charge
02807     // scheme is being used.
02808     QMAtomData *pcPpme = NULL;
02809     if (msg->switching) {
02810         
02811         if (msg->PMEOn)
02812             pcPpme = new QMAtomData[msg->numRealPntChrgs];
02813         
02814         pntChrgSwitching(msg, pcPpme) ;
02815     } else {
02816         pcPpme = pcP;
02817     }
02818     
02819     int retVal = 0;
02820     struct stat info;
02821     
02822     // For each QM group, create a subdirectory where files will be palced.
02823     std::string baseDir(msg->baseDir);
02824     std::ostringstream itosConv ;
02825     if ( CmiNumPartitions() > 1 ) {
02826         baseDir.append("/") ;
02827         itosConv << CmiMyPartition() ;
02828         baseDir += itosConv.str() ;
02829         itosConv.str("");
02830         itosConv.clear() ;
02831         
02832         if (stat(msg->baseDir, &info) != 0 ) {
02833             CkPrintf( "Node %d cannot access directory %s\n",
02834                       CkMyPe(), baseDir.c_str() );
02835             NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
02836         }
02837         else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
02838             DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
02839             retVal = mkdir(baseDir.c_str(), S_IRWXU);
02840         }
02841         
02842     }
02843     baseDir.append("/") ;
02844     itosConv << msg->peIter ;
02845     baseDir += itosConv.str() ;
02846     
02847     if (stat(msg->baseDir, &info) != 0 ) {
02848         CkPrintf( "Node %d cannot access directory %s\n",
02849                   CkMyPe(), baseDir.c_str() );
02850         NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
02851     }
02852     else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
02853         DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
02854         retVal = mkdir(baseDir.c_str(), S_IRWXU);
02855     }
02856     
02857     #ifdef DEBUG_QM
02858     Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
02859     #endif
02860     
02861     inputFileName.clear();
02862     inputFileName.append(baseDir.c_str()) ;
02863     inputFileName.append("/qmmm_") ;
02864     inputFileName += itosConv.str() ;
02865     inputFileName.append(".input") ;
02866     
02867     // Opens file for coordinate and parameter input
02868     inputFile = fopen(inputFileName.c_str(),"w");
02869     if ( ! inputFile ) {
02870         iout << iERROR << "Could not open input file for writing: " 
02871         << inputFileName << "\n" << endi ;
02872         NAMD_err("Error writing QM input file.");
02873     }
02874     
02875     // Builds the command that will be executed
02876     qmCommand.clear();
02877     qmCommand.append("cd ");
02878     qmCommand.append(baseDir);
02879     qmCommand.append(" ; ");
02880     qmCommand.append(msg->execPath) ;
02881     qmCommand.append(" ") ;
02882     qmCommand.append(inputFileName) ;
02883     qmCommand.append(" > /dev/null 2>&1") ;
02884     
02885     // Builds the file name where MOPAC will place the gradient
02886     // This is also relative to the input file name.
02887     outputFileName = inputFileName ;
02888     outputFileName.append(".aux") ;
02889     
02890     if (msg->numAllPntChrgs) {
02891         // Builds the file name where we will write the point charges.
02892         pntChrgFileName.clear();
02893         pntChrgFileName.append(baseDir.c_str()) ;
02894         pntChrgFileName.append("/mol.in") ;
02895         
02896         chrgFile = fopen(pntChrgFileName.c_str(),"w");
02897         if ( ! chrgFile ) {
02898             iout << iERROR << "Could not open charge file for writing: " 
02899             << pntChrgFileName << "\n" << endi ;
02900             NAMD_err("Error writing point charge file.");
02901         }
02902         
02903         iret = fprintf(chrgFile,"\n%d %d\n", 
02904                        msg->numQMAtoms, msg->numAllAtoms - msg->numQMAtoms);
02905         if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
02906     }
02907     
02908     // writes configuration lines to the input file.
02909     std::stringstream ss(msg->configLines) ;
02910     std::string confLineStr;
02911     while (std::getline(ss, confLineStr) ) {
02912         confLineStr.append("\n");
02913         iret = fprintf(inputFile,confLineStr.c_str());
02914         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
02915     }
02916     
02917     
02918     iret = fprintf(inputFile,"\n");
02919     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
02920     
02921     DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " 
02922     << inputFileName.c_str() << " and " << msg->numAllPntChrgs 
02923     << " point charges in file " << pntChrgFileName.c_str() << "\n");
02924     
02925     // write QM and dummy atom coordinates to input file and
02926     // MM electric field from MM point charges.
02927     for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
02928         
02929         double x = atmP->position.x;
02930         double y = atmP->position.y;
02931         double z = atmP->position.z;
02932         
02933         iret = fprintf(inputFile,"%s %f 1 %f 1 %f 1\n",
02934                        atmP->element,x,y,z);
02935         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
02936         
02937         if (msg->numAllPntChrgs) {
02938             BigReal phi = 0;
02939             
02940             // The Electrostatic Potential is calculated according to
02941             // the QMMM Keyword documentation for MOPAC
02942             // http://openmopac.net/manual/QMMM.html
02943             pcP = msg->data + msg->numAllAtoms ;
02944             for ( size_t j=0; j < msg->numAllPntChrgs; ++j, ++pcP ) {
02945                 
02946                 // In case of QM-MM bonds, the charge of the MM1 atom is ignored
02947                 if (pcP->type == QMPCTYPE_IGNORE)
02948                     continue;
02949                 
02950                 double charge = pcP->charge;
02951                 
02952                 double xMM = pcP->position.x;
02953                 double yMM = pcP->position.y;
02954                 double zMM = pcP->position.z;
02955                 
02956                 BigReal rij = sqrt((x-xMM)*(x-xMM) +
02957                      (y-yMM)*(y-yMM) +
02958                      (z-zMM)*(z-zMM) ) ;
02959                 
02960                 phi += charge/rij ;
02961             }
02962             
02963             // We use the same Coulomb constant used in the rest of NAMD
02964             // instead of the suggested rounded 332 suggested by MOPAC.
02965             phi = phi*constants ;
02966             
02967             iret = fprintf(chrgFile,"%s %f %f %f %f\n",
02968                            atmP->element,x,y,z, phi);
02969             if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
02970         }
02971     }
02972     
02973     DebugM(4,"Closing input and charge file\n");
02974     
02975     if (msg->numAllPntChrgs)
02976         fclose(chrgFile);
02977     
02978     fclose(inputFile);
02979     
02980     if (msg->prepProcOn) {
02981         
02982         std::string prepProc(msg->prepProc) ;
02983         prepProc.append(" ") ;
02984         prepProc.append(inputFileName) ;
02985         iret = system(prepProc.c_str());
02986         if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
02987         if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
02988     }
02989     
02990     // runs QM command
02991     DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
02992     iret = system(qmCommand.c_str());
02993     
02994     if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
02995     if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
02996 
02997     if (msg->secProcOn) {
02998         
02999         std::string secProc(msg->secProc) ;
03000         secProc.append(" ") ;
03001         secProc.append(inputFileName) ;
03002         itosConv.str("");
03003         itosConv.clear() ;
03004         itosConv << msg->timestep ;
03005         secProc.append(" ") ;
03006         secProc += itosConv.str() ;
03007         
03008         iret = system(secProc.c_str());
03009         if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
03010         if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
03011     }
03012     
03013     // remove coordinate file
03014 //     iret = remove(inputFileName);
03015 //     if ( iret ) { NAMD_err(0); }
03016 
03017     // remove coordinate file
03018 //     iret = remove(pntChrgFileName);
03019 //     if ( iret ) { NAMD_err(0); }
03020     
03021     // opens output file
03022     DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
03023     outputFile = fopen(outputFileName.c_str(),"r");
03024     if ( ! outputFile ) {
03025         iout << iERROR << "Could not find QM output file!\n" << endi;
03026         NAMD_err(0); 
03027     }
03028     
03029     // Resets the pointers.
03030     atmP = msg->data ;
03031     pcP = msg->data + msg->numAllAtoms ;
03032     
03033     // Allocates the return message.
03034     QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
03035     resMsg->grpIndx = msg->grpIndx;
03036     resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
03037     resMsg->energyOrig = 0;
03038     resMsg->energyCorr = 0;
03039     for ( int k=0; k<3; ++k )
03040         for ( int l=0; l<3; ++l )
03041             resMsg->virial[k][l] = 0;
03042     QMForce *resForce = resMsg->force ;
03043     
03044     // We split the data into two pointers so we can "skip" the dummy atoms
03045     // which have no representation in NAMD.
03046     for (int i=0; i<resMsg->numForces; i++) {
03047         resForce[i].force = 0;
03048         resForce[i].charge = 0 ;
03049         if (i < msg->numQMAtoms) {
03050             // We use the replace field to indicate QM atoms,
03051             // which will have charge information.
03052             resForce[i].replace = 1 ;
03053             resForce[i].id = atmP->id;
03054             atmP++;
03055         }
03056         else {
03057             // We use the replace field to indicate QM atoms,
03058             // which will have charge information.
03059             resForce[i].replace = 0 ;
03060             resForce[i].id = pcP->id;
03061             pcP++;
03062         }
03063     }
03064     
03065     // Resets the pointers.
03066     atmP = msg->data ;
03067     pcP = msg->data + msg->numAllAtoms ;
03068     
03069     // Reads the data form the output file created by the QM software.
03070     // Gradients over the QM atoms, and Charges for QM atoms will be read.
03071     
03072     size_t atmIndx = 0, gradCount = 0;
03073     Bool gradFields = false, chargeFields = false;
03074     Bool chargesRead = false, gradsRead = false;
03075     while ( fgets(line, lineLen, outputFile) != NULL){
03076         // We loop over the lines of the output file untill we find
03077         // the line that initiates the "atom charges" lines. Then
03078         // we read all lines and expect to get one or more charges 
03079         // per line, separated by space(s), untill we find a line that
03080         // initiates the "gradients", then once more, we expect several
03081         // numbers separated by space(s). When the "overlap matrix" 
03082         // string is found, we break the loop and stop reading the file.
03083         
03084         // if ( strstr(line,"TOTAL_ENERGY") != NULL ) {
03085         if ( strstr(line,"HEAT_OF_FORMATION") != NULL ) {
03086             
03087             char strEnergy[14], *endPtr ;
03088             
03089             //strncpy(strEnergy, line + 17, 13) ;
03090             strncpy(strEnergy, line + 28, 13) ;
03091             strEnergy[13] = '\0';
03092             
03093             // We have to convert the notation from FORTRAN double precision to
03094             // the natural, normal, reasonable and not terribly out dated format.
03095             resMsg->energyOrig = strtod(strEnergy, &endPtr);
03096             if ( *endPtr == 'D' ) {
03097                 *endPtr = 'e';
03098                 resMsg->energyOrig = strtod(strEnergy, &endPtr);
03099             }
03100             
03101             // In MOPAC, the total energy is given in EV, so we convert to Kcal/mol
03102 //             resMsg->energyOrig *= 23.060 ; // We now read Heat of Formation, which is given in Kcal/Mol
03103             
03104 //             DebugM(4,"Reading QM energy from file: " << resMsg->energyOrig << "\n");
03105             
03106             resMsg->energyCorr = resMsg->energyOrig;
03107             
03108             continue;
03109         }
03110         
03111         if ( strstr(line,"ATOM_CHARGES") != NULL ) {
03112             gradFields = false;
03113             chargeFields = true;
03114             atmIndx = 0;
03115             
03116             // If the user asked for charges NOT to be read from MOPAC,
03117             // we skip the charge reading step.
03118             if (msg->qmAtmChrgMode == QMCHRGNONE) {
03119                 chargeFields = false;
03120                 atmIndx = msg->numAllAtoms;
03121             }
03122             else {
03123                 chargeFields = true;
03124                 atmIndx = 0;
03125             }
03126             
03127             // Now we expect the following line(s) to have atom charges
03128             continue;
03129         }
03130         
03131         if ( strstr(line,"GRADIENTS") != NULL ) {
03132             
03133             // Now that all charges have been read, checks if the 
03134             // numbers match
03135             if (atmIndx != msg->numAllAtoms) {
03136                 NAMD_die("Error reading QM forces file. Wrong number of atom charges");
03137             }
03138             
03139             chargesRead = true;
03140             
03141             // Now we expect the following line(s) to have gradients
03142             chargeFields = false ;
03143             gradFields = true;
03144             gradCount = 0;
03145             atmIndx = 0;
03146             
03147             continue;
03148         }
03149         
03150         if ( strstr(line,"OVERLAP_MATRIX") != NULL ) {
03151             
03152             // Now that all gradients have been read, checks if the 
03153             // numbers match
03154             if (atmIndx != msg->numAllAtoms) {
03155                 NAMD_die("Error reading QM forces file. Wrong number of gradients");
03156             }
03157             
03158             gradsRead = true;
03159             
03160             // Nothing more to read from the ".aux" file
03161             break;
03162         }
03163         
03164         char result[10] ;
03165         size_t strIndx = 0;
03166         
03167         if (chargeFields) {
03168             while ((strIndx < (strlen(line)-9)) && (strlen(line)-1 >=9 ) ) {
03169                 
03170                 strncpy(result, line+strIndx,9) ;
03171                 result[9] = '\0';
03172                 
03173                 Real localCharge = atof(result);
03174                 
03175                 // If we are reading charges from QM atoms, store them.
03176                 if (atmIndx < msg->numQMAtoms ) {
03177                     atmP[atmIndx].charge = localCharge;
03178                     resForce[atmIndx].charge = localCharge;
03179                 }
03180                 
03181                 // If we are reading charges from Dummy atoms,
03182                 // place them on the appropriate QM atoms.
03183                 if ( msg->numQMAtoms <= atmIndx &&
03184                     atmIndx < msg->numAllAtoms ) {
03185                     // We keep the calculated charge in the dummy atom
03186                     // so that we can calculate electrostatic forces later on.
03187                     atmP[atmIndx].charge = localCharge;
03188                     
03189                     // The dummy atom points to the QM atom to which it is bound.
03190                     int qmInd = atmP[atmIndx].bountToIndx ;
03191                     resForce[qmInd].charge += localCharge;
03192                 }
03193                 
03194                 strIndx += 9;
03195                 atmIndx++;
03196                 
03197                 // If we found all charges for QM and dummy atoms, break the loop
03198                 // and stop reading the line.
03199                 if (atmIndx == msg->numAllAtoms) {
03200                     chargeFields = false;
03201                     break;
03202                 }
03203                 
03204             }
03205             
03206         }
03207         
03208         if (gradFields) {
03209             while ((strIndx < (strlen(line)-9)) && (strlen(line)-1 >=9 ) ) {
03210                 
03211                 strncpy(result, line+strIndx,9) ;
03212                 result[9] = '\0';
03213                 
03214                 gradient[gradCount] = atof(result);
03215                 if (gradCount == 2) {
03216                     
03217                     if (atmIndx < msg->numQMAtoms ) {
03218                         // Gradients in MOPAC are written in kcal/mol/A.
03219                         resForce[atmIndx].force.x = -1*gradient[0];
03220                         resForce[atmIndx].force.y = -1*gradient[1];
03221                         resForce[atmIndx].force.z = -1*gradient[2];
03222                     }
03223                     
03224                     // If we are reading forces applied on Dummy atoms,
03225                     // place them on the appropriate QM and MM atoms to conserve energy.
03226                     
03227                     // This implementation was based on the description in 
03228                     // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
03229                     if ( msg->numQMAtoms <= atmIndx &&
03230                     atmIndx < msg->numAllAtoms ) {
03231                         // The dummy atom points to the QM atom to which it is bound.
03232                         // The QM atom and the MM atom (in a QM-MM bond) point to each other.
03233                         int qmInd = atmP[atmIndx].bountToIndx ;
03234                         int mmInd = atmP[qmInd].bountToIndx ;
03235                         
03236                         Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
03237                         Real mmqmDist = dir.length() ;
03238                         
03239                         Real linkDist = Vector(atmP[atmIndx].position - 
03240                                         atmP[qmInd].position).length() ;
03241                         
03242                         Force mmForce(0), qmForce(0), 
03243                             linkForce(gradient[0], gradient[1], gradient[2]);
03244                         linkForce *= -1;
03245                         
03246                         Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
03247                         // Unit vectors
03248                         Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
03249                         Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
03250                         Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
03251                         Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
03252                         
03253                         qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv + 
03254                                     (xDelta)*base) )*xuv;
03255                         
03256                         qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv + 
03257                                     (yDelta)*base) )*yuv;
03258                         
03259                         qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv + 
03260                                     (zDelta)*base) )*zuv;
03261                         
03262                         
03263                         mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
03264                                     (xDelta)*base) )*xuv;
03265                         
03266                         mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
03267                                     (yDelta)*base) )*yuv;
03268                         
03269                         mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
03270                                     (zDelta)*base) )*zuv;
03271                         
03272                         resForce[qmInd].force += qmForce;
03273                         resForce[msg->numQMAtoms + mmInd].force += mmForce;
03274                     }
03275                     
03276                     gradCount = 0;
03277                     atmIndx++;
03278                 } else {
03279                     gradCount++;
03280                 }
03281                 
03282                 strIndx += 9;
03283                 
03284                 // If we found all gradients for QM atoms, break the loop
03285                 // and keep the next gradient line from being read, if there
03286                 // is one. Following gradients, if present, will refer to link
03287                 // atoms, and who cares about those?.
03288                 if (atmIndx == msg->numAllAtoms) {
03289                     gradFields = false;
03290                     break;
03291                 }
03292             }
03293         }
03294         
03295     }
03296     
03297     delete [] line;
03298     
03299     fclose(outputFile);
03300     
03301     // In case charges are not to be read form the QM software output,
03302     // we load the origianl atom charges.
03303     if (msg->qmAtmChrgMode == QMCHRGNONE) {
03304         int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
03305         const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
03306         Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
03307         
03308         atmIndx = 0 ;
03309         for (; atmIndx < msg->numQMAtoms; atmIndx++) {
03310             
03311             // Loops over all QM atoms (in all QM groups) comparing their global indices
03312             for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
03313                 
03314                 if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
03315                     
03316                     atmP[atmIndx].charge = qmAtmChrg[qmIter];
03317                     resForce[atmIndx].charge = qmAtmChrg[qmIter];
03318                     
03319                     break;
03320                 }
03321                 
03322             }
03323             
03324         }
03325         for (size_t j=msg->numQMAtoms; j<msg->numAllAtoms; ++j ) {
03326             atmP[j].charge = 0;
03327         }
03328     }
03329     
03330     // remove force file
03331 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
03332 //     iret = remove(outputFileName);
03333 //     if ( iret ) { NAMD_err(0); }
03334     
03335     if (! (chargesRead && gradsRead) ) {
03336         NAMD_die("Error reading QM forces file. Not all data could be read!");
03337     }
03338     
03339     DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
03340     
03341     atmP = msg->data ;
03342     pcP = msg->data + msg->numAllAtoms ;
03343     
03344     // The initial point charge index for the *force* message is the number of QM
03345     // atoms, since the dummy atoms have no representation in NAMD
03346     int pcIndx = msg->numQMAtoms;
03347     
03348     // ---- WARNING  ----
03349             // We need to apply the force felt by each QM atom due to the classical 
03350             // point charges sent to MOPAC.
03351             // MOPAC takes the MM electrostatic potential into cosideration ONLY
03352             // when performing the SCF calculation and when returning the energy
03353             // of the system, but it does not apply the potential to the final
03354             // gradient calculation, therefore, we calculate the resulting force
03355             // here.
03356     // Initialize force vector for electrostatic contribution
03357     std::vector<Force> qmElecForce ;
03358     for (size_t j=0; j<msg->numAllAtoms; ++j )
03359         qmElecForce.push_back( Force(0) );
03360     
03361     // We loop over point charges and distribute the forces applied over
03362     // virtual point charges to the MM1 and MM2 atoms (if any virtual PCs exists).
03363     for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
03364         
03365         // No force was applied to the QM region due to this charge, so it 
03366         // does not receive any back from the QM region. It must be an MM1 atom
03367         // from a QM-MM bond.
03368         if (pcP[i].type == QMPCTYPE_IGNORE)
03369             continue;
03370         
03371         Force totalForce(0);
03372         
03373         BigReal pntCharge = pcP[i].charge;
03374         
03375         Position posMM = pcP[i].position ;
03376         
03377         for (size_t j=0; j<msg->numAllAtoms; ++j ) {
03378             
03379             BigReal qmCharge = atmP[j].charge ;
03380             
03381             BigReal force = pntCharge*qmCharge*constants ;
03382             
03383             Position rVec = posMM - atmP[j].position ;
03384             
03385             force /= rVec.length2();
03386             
03387             // We accumulate the total force felt by a point charge
03388             // due to the QM charge distribution. This total force
03389             // will be applied to the point charge if it is a real one,
03390             // or will be distirbuted over MM1 and MM2 point charges, it 
03391             // this is a virtual point charge.
03392             totalForce += force*rVec.unit();
03393             
03394             // Accumulate force felt by a QM atom due to point charges.
03395             qmElecForce[j] += -1*force*rVec.unit();
03396         }
03397         
03398         if (pcP[i].type == QMPCTYPE_CLASSICAL) {
03399             // If we already ignored MM1 charges, we take all other 
03400             // non-virtual charges and apply forces directly to them.
03401             resForce[pcIndx].force += totalForce;
03402         }
03403         else {
03404             // If we are handling virtual PC, we distribute the force over
03405             // MM1 and MM2.
03406             
03407             // Virtual PC are bound to MM2.
03408             int mm2Indx = pcP[i].bountToIndx;
03409             // MM2 charges are bound to MM1.
03410             int mm1Indx = pcP[mm2Indx].bountToIndx;
03411             
03412             Real Cq = pcP[i].dist;
03413             
03414             Force mm1Force = (1-Cq)*totalForce ;
03415             Force mm2Force = Cq*totalForce ;
03416             
03417             resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
03418             resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
03419         }
03420         
03421     }
03422     
03423     // We now apply the accumulated electrostatic force contribution
03424     // to QM atoms.
03425     for (size_t j=0; j<msg->numAllAtoms; ++j ) {
03426         
03427         if (j < msg->numQMAtoms ) {
03428             resForce[j].force += qmElecForce[j];
03429             
03430         } else {
03431             // In case the QM atom is a Link atom, we redistribute 
03432             // its force as bove.
03433             
03434             // The dummy atom points to the QM atom to which it is bound.
03435             // The QM atom and the MM atom (in a QM-MM bond) point to each other.
03436             int qmInd = atmP[j].bountToIndx ;
03437             int mmInd = atmP[qmInd].bountToIndx ;
03438             
03439             Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
03440             Real mmqmDist = dir.length() ;
03441             
03442             Real linkDist = Vector(atmP[j].position - 
03443                             atmP[qmInd].position).length() ;
03444             
03445             Force linkForce = qmElecForce[j];
03446             
03447             Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
03448             // Unit vectors
03449             Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
03450             Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
03451             Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
03452             Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
03453             
03454             Force qmForce = (linkForce*((1 - linkDist/mmqmDist)*xuv + 
03455                         (xDelta)*base) )*xuv;
03456             
03457             qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv + 
03458                         (yDelta)*base) )*yuv;
03459             
03460             qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv + 
03461                         (zDelta)*base) )*zuv;
03462             
03463             
03464             Force mmForce = (linkForce*((linkDist/mmqmDist)*xuv -
03465                         (xDelta)*base) )*xuv;
03466             
03467             mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
03468                         (yDelta)*base) )*yuv;
03469             
03470             mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
03471                         (zDelta)*base) )*zuv;
03472             
03473             resForce[qmInd].force += qmForce;
03474             resForce[msg->numQMAtoms + mmInd].force += mmForce;
03475         }
03476     }
03477     
03478     // Adjusts forces from PME, canceling contributions from the QM and 
03479     // direct Coulomb forces calculated here.
03480     if (msg->PMEOn) {
03481         
03482         DebugM(1,"Correcting forces and energy for PME.\n");
03483         
03484         ewaldcof = msg->PMEEwaldCoefficient;
03485         BigReal TwoBySqrtPi = 1.12837916709551;
03486         pi_ewaldcof = TwoBySqrtPi * ewaldcof;
03487         
03488         for (size_t i=0; i < msg->numQMAtoms; i++) {
03489             
03490             BigReal p_i_charge = atmP[i].charge ;
03491             Position pos_i = atmP[i].position ;
03492             
03493             const BigReal kq_i = p_i_charge * constants;
03494             
03495             for (size_t j=i+1; j < msg->numQMAtoms; j++) {
03496                 
03497                 BigReal p_j_charge = atmP[j].charge ;
03498                 
03499                 Position pos_j = atmP[j].position ;
03500                 
03501                 BigReal r = Vector(pos_i - pos_j).length() ;
03502                 
03503                 BigReal tmp_a = r * ewaldcof;
03504                 BigReal tmp_b = erfc(tmp_a);
03505                 BigReal corr_energy = tmp_b;
03506                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
03507                 
03508 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
03509                 BigReal recip_energy = (1-tmp_b)/r;
03510                 
03511 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
03512                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
03513                 
03514                 // Final force and energy correction for this pair of atoms.
03515                 BigReal energy = kq_i * p_j_charge * recip_energy ;
03516                 
03517                 Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
03518                 
03519                 // The force is *subtracted* from the total force acting on
03520                 // both atoms. The sign on fixForce corrects the orientation
03521                 // of the subtracted force.
03522 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
03523 //                     << std::endl);
03524 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
03525 //                     << std::endl);
03526 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
03527 //                 DebugM(4,"Energy correction: " << energy << "\n");
03528                 resForce[i].force -= fixForce ;
03529                 resForce[j].force -= -1*fixForce ;
03530                 
03531                 // The energy is *subtracted* from the total energy calculated here.
03532                 resMsg->energyCorr -= energy;
03533             }
03534         }
03535         
03536 //         DebugM(4,"Corrected energy QM-QM interactions: " << resMsg->energyCorr << "\n");
03537         
03538         pcIndx = msg->numQMAtoms;
03539         // We only loop over point charges from real atoms, ignoring the ones 
03540         // created to handle QM-MM bonds.
03541         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
03542             
03543             BigReal p_i_charge = pcPpme[i].charge;
03544             Position pos_i = pcPpme[i].position ;
03545             
03546             const BigReal kq_i = p_i_charge * constants;
03547             
03548             Force fixForce = 0;
03549             
03550             for (size_t j=0; j<msg->numQMAtoms; ++j ) {
03551                 
03552                 BigReal p_j_charge = atmP[j].charge ;
03553                 
03554                 Position pos_j = atmP[j].position ;
03555                 
03556                 BigReal r = Vector(pos_i - pos_j).length() ;
03557                 
03558                 BigReal tmp_a = r * ewaldcof;
03559                 BigReal tmp_b = erfc(tmp_a);
03560                 BigReal corr_energy = tmp_b;
03561                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
03562                 
03563 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
03564                 BigReal recip_energy = (1-tmp_b)/r;
03565                 
03566 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
03567                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
03568                 
03569                 // Final force and energy correction for this pair of atoms.
03570                 BigReal energy = kq_i * p_j_charge * recip_energy ;
03571                 
03572                 fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
03573                 
03574                 // The energy is *subtracted* from the total energy calculated here.
03575                 resMsg->energyCorr -= energy;
03576             }
03577             
03578             // The force is *subtracted* from the total force acting on
03579             // the point charge..
03580 //             DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
03581 //             << std::endl);
03582 //             DebugM(4,"Force correction: " << fixForce << std::endl);
03583             resForce[pcIndx].force -= kq_i*fixForce ;
03584         }
03585         
03586     }
03587     
03588     // Calculates the virial contribution form the forces on QM atoms and 
03589     // point charges calculated here.
03590     for (size_t i=0; i < msg->numQMAtoms; i++) {
03591         // virial used by NAMD is -'ve of normal convention, so reverse it!
03592         // virial[i][j] in file should be sum of -1 * f_i * r_j
03593         for ( int k=0; k<3; ++k )
03594             for ( int l=0; l<3; ++l )
03595                 resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
03596     }
03597     
03598     pcIndx = msg->numQMAtoms; // Index in the real PC force array.
03599     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
03600         // virial used by NAMD is -'ve of normal convention, so reverse it!
03601         // virial[i][j] in file should be sum of -1 * f_i * r_j
03602         for ( int k=0; k<3; ++k )
03603             for ( int l=0; l<3; ++l )
03604                 resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
03605     }
03606     
03607     
03608     
03609     // Send message to rank zero with results.
03610     QMProxy[0].recvQMRes(resMsg);
03611     
03612     if (msg->switching && msg->PMEOn)
03613         delete [] pcPpme;
03614     
03615     delete msg;
03616     return ;
03617 }

void ComputeQMMgr::calcORCA ( QMGrpCalcMsg  ) 

Definition at line 3619 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMForce::charge, QMAtomData::charge, QMGrpCalcMsg::charge, QMGrpCalcMsg::configLines, QMGrpCalcMsg::constants, QMGrpCalcMsg::data, DebugM, endi(), QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, QMGrpCalcMsg::execPath, QMGrpResMsg::force, QMGrpCalcMsg::grpIndx, QMGrpResMsg::grpIndx, QMAtomData::id, iERROR(), if(), iout, j, Vector::length(), QMGrpCalcMsg::multiplicity, NAMD_die(), NAMD_err(), QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMGrpResMsg::numForces, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, QMGrpCalcMsg::peIter, QMGrpCalcMsg::PMEOn, QMAtomData::position, QMGrpCalcMsg::prepProc, QMGrpCalcMsg::prepProcOn, QMPCTYPE_IGNORE, QMGrpCalcMsg::secProc, QMGrpCalcMsg::secProcOn, QMGrpCalcMsg::switching, QMGrpCalcMsg::timestep, QMAtomData::type, QMGrpResMsg::virial, Vector::x, x, Vector::y, y, Vector::z, and z.

03620 {
03621     
03622     FILE *inputFile,*outputFile,*chrgFile;
03623     int iret;
03624     
03625     const size_t lineLen = 256;
03626     char *line = new char[lineLen];
03627     
03628     std::string qmCommand, inputFileName, outputFileName, pntChrgFileName, pcGradFilename;
03629     std::string tmpRedirectFileName, pcGradFileName;
03630     
03631     // For coulomb calculation
03632     BigReal constants = msg->constants;
03633     
03634     double gradient[3];
03635     
03636     DebugM(4,"Running ORCA on PE " << CkMyPe() << std::endl);
03637     
03638     QMAtomData *pcP = msg->data + msg->numAllAtoms ;
03639     
03640     // We create a pointer for PME correction, which may point to
03641     // a copy of the original point charge array, with unchanged charges,
03642     // or points to the original array in case no switching or charge
03643     // scheme is being used.
03644     QMAtomData *pcPpme = NULL;
03645     if (msg->switching) {
03646         
03647         if (msg->PMEOn)
03648             pcPpme = new QMAtomData[msg->numRealPntChrgs];
03649         
03650         pntChrgSwitching(msg, pcPpme) ;
03651     } else {
03652         pcPpme = pcP;
03653     }
03654     
03655     // For each QM group, create a subdirectory where files will be palced.
03656     std::string baseDir(msg->baseDir);
03657     baseDir.append("/") ;
03658     std::ostringstream itosConv ;
03659     itosConv << msg->peIter ;
03660     baseDir += itosConv.str() ;
03661     
03662     struct stat info;
03663     
03664     if (stat(msg->baseDir, &info) != 0 ) {
03665         CkPrintf( "Node %d cannot access directory %s\n",
03666                   CkMyPe(), baseDir.c_str() );
03667         NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
03668     }
03669     else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
03670         DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
03671         int retVal = mkdir(baseDir.c_str(), S_IRWXU);
03672     }
03673     
03674     #ifdef DEBUG_QM
03675     Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
03676     #endif
03677     
03678     inputFileName.clear();
03679     inputFileName.append(baseDir.c_str()) ;
03680     inputFileName.append("/qmmm_") ;
03681     inputFileName += itosConv.str() ;
03682     inputFileName.append(".input") ;
03683     
03684     // Opens file for coordinate and parameter input
03685     inputFile = fopen(inputFileName.c_str(),"w");
03686     if ( ! inputFile ) {
03687         iout << iERROR << "Could not open input file for writing: " 
03688         << inputFileName << "\n" << endi ;
03689         NAMD_err("Error writing QM input file.");
03690     }
03691     
03692     // Builds the command that will be executed
03693     qmCommand.clear();
03694     qmCommand.append("cd ");
03695     qmCommand.append(baseDir);
03696     qmCommand.append(" ; ");
03697     qmCommand.append(msg->execPath) ;
03698     qmCommand.append(" ") ;
03699     qmCommand.append(inputFileName) ;
03700     
03701     // Build the redirect file name we need for the screen output.
03702     // That's the only place where we can get partial charges for QM atoms.
03703     tmpRedirectFileName = inputFileName ;
03704     tmpRedirectFileName.append(".TmpOut") ;
03705     
03706     qmCommand.append(" > ") ;
03707     qmCommand.append(tmpRedirectFileName) ;
03708     
03709     // Builds the file name where orca will place the gradient
03710     // This will be relative to the input file
03711     outputFileName = inputFileName ;
03712     outputFileName.append(".engrad") ;
03713     
03714     // Builds the file name where orca will place gradients acting on 
03715     // point charges
03716     pcGradFilename = inputFileName ;
03717     pcGradFilename.append(".pcgrad") ;
03718     
03719     if (msg->numAllPntChrgs) {
03720         // Builds the file name where we will write the point charges.
03721         pntChrgFileName = inputFileName ;
03722         pntChrgFileName.append(".pntchrg") ;
03723         
03724         pcGradFileName = inputFileName;
03725         pcGradFileName.append(".pcgrad");
03726         
03727         chrgFile = fopen(pntChrgFileName.c_str(),"w");
03728         if ( ! chrgFile ) {
03729             iout << iERROR << "Could not open charge file for writing: " 
03730             << pntChrgFileName << "\n" << endi ;
03731             NAMD_err("Error writing point charge file.");
03732         }
03733         
03734         int numPntChrgs = 0;
03735         for (int i=0; i<msg->numAllPntChrgs; i++ ) {
03736             if (pcP[i].type != QMPCTYPE_IGNORE)
03737                 numPntChrgs++;
03738         }
03739         
03740         iret = fprintf(chrgFile,"%d\n", numPntChrgs);
03741         if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
03742     }
03743     
03744     // writes configuration lines to the input file.
03745     std::stringstream ss(msg->configLines) ;
03746     std::string confLineStr;
03747     while (std::getline(ss, confLineStr) ) {
03748         confLineStr.append("\n");
03749         iret = fprintf(inputFile,confLineStr.c_str());
03750         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03751     }
03752     
03753     if (msg->numAllPntChrgs) {
03754         iret = fprintf(inputFile,"%%pointcharges \"%s\"\n", pntChrgFileName.c_str());
03755         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03756     }
03757     
03758     iret = fprintf(inputFile,"\n\n%%coords\n  CTyp xyz\n");
03759     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03760     
03761     iret = fprintf(inputFile,"  Charge %f\n",msg->charge);
03762     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03763     
03764     iret = fprintf(inputFile,"  Mult %f\n",msg->multiplicity);
03765     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03766     
03767     iret = fprintf(inputFile,"  Units Angs\n  coords\n\n");
03768     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03769     
03770     DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " << 
03771     inputFileName.c_str() << " and " << msg->numAllPntChrgs << 
03772     " point charges in file " << pntChrgFileName.c_str() << "\n");
03773     
03774     // write QM and dummy atom coordinates to input file.
03775     QMAtomData *atmP = msg->data ;
03776     for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
03777         
03778         double x = atmP->position.x;
03779         double y = atmP->position.y;
03780         double z = atmP->position.z;
03781         
03782         iret = fprintf(inputFile,"  %s %f %f %f\n",
03783                        atmP->element,x,y,z);
03784         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03785         
03786     }
03787     
03788     iret = fprintf(inputFile,"  end\nend\n");
03789     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
03790     
03791     if (msg->numAllPntChrgs) {
03792         // Write point charges to file.
03793         pcP = msg->data + msg->numAllAtoms ;
03794         for ( size_t j=0; j < msg->numAllPntChrgs; j++, ++pcP) {
03795             
03796             if (pcP->type == QMPCTYPE_IGNORE)
03797                 continue;
03798             
03799             double charge = pcP->charge;
03800             
03801             double x = pcP->position.x;
03802             double y = pcP->position.y;
03803             double z = pcP->position.z;
03804             
03805             iret = fprintf(chrgFile,"%f %f %f %f\n",
03806                            charge,x,y,z);
03807             if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
03808         }
03809         
03810         fclose(chrgFile);
03811     }
03812     
03813     DebugM(4,"Closing input and charge file\n");
03814     fclose(inputFile);
03815     
03816     if (msg->prepProcOn) {
03817         
03818         std::string prepProc(msg->prepProc) ;
03819         prepProc.append(" ") ;
03820         prepProc.append(inputFileName) ;
03821         iret = system(prepProc.c_str());
03822         if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
03823         if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
03824     }
03825     
03826         // runs QM command
03827     DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
03828     iret = system(qmCommand.c_str());
03829     
03830     if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
03831     if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
03832 
03833     if (msg->secProcOn) {
03834         
03835         std::string secProc(msg->secProc) ;
03836         secProc.append(" ") ;
03837         secProc.append(inputFileName) ;
03838         itosConv.str("");
03839         itosConv.clear() ;
03840         itosConv << msg->timestep ;
03841         secProc.append(" ") ;
03842         secProc += itosConv.str() ;
03843         
03844         iret = system(secProc.c_str());
03845         if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
03846         if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
03847     }
03848 
03849     // remove coordinate file
03850 //     iret = remove(inputFileName);
03851 //     if ( iret ) { NAMD_err(0); }
03852 
03853     // remove coordinate file
03854 //     iret = remove(pntChrgFileName);
03855 //     if ( iret ) { NAMD_err(0); }
03856     
03857     // opens output file
03858     DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
03859     outputFile = fopen(outputFileName.c_str(),"r");
03860     if ( ! outputFile ) {
03861         iout << iERROR << "Could not find QM output file!\n" << endi;
03862         NAMD_err(0); 
03863     }
03864 
03865     // Resets the pointers.
03866     atmP = msg->data ;
03867     pcP = msg->data + msg->numAllAtoms ;
03868     
03869     // Allocates the return message.
03870     QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
03871     resMsg->grpIndx = msg->grpIndx;
03872     resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
03873     resMsg->energyOrig = 0;
03874     resMsg->energyCorr = 0;
03875     for ( int k=0; k<3; ++k )
03876         for ( int l=0; l<3; ++l )
03877             resMsg->virial[k][l] = 0;
03878     QMForce *resForce = resMsg->force ;
03879     
03880     // We split the data into two pointers so we can "skip" the dummy atoms
03881     // which have no representation in NAMD.
03882     for (int i=0; i<resMsg->numForces; i++) {
03883         resForce[i].force = 0;
03884         resForce[i].charge = 0 ;
03885         if (i < msg->numQMAtoms) {
03886             // We use the replace field to indicate QM atoms,
03887             // which will have charge information.
03888             resForce[i].replace = 1 ;
03889             resForce[i].id = atmP->id;
03890             atmP++;
03891         }
03892         else {
03893             // We use the replace field to indicate QM atoms,
03894             // which will have charge information.
03895             resForce[i].replace = 0 ;
03896             resForce[i].id = pcP->id;
03897             pcP++;
03898         }
03899     }
03900     
03901     // Resets the pointers.
03902     atmP = msg->data ;
03903     pcP = msg->data + msg->numAllAtoms ;
03904     
03905     size_t atmIndx = 0, gradCount = 0;
03906     int numAtomsInFile = 0;
03907     
03908     // Reads the data form the output file created by the QM software.
03909     // Gradients over the QM atoms, and Charges for QM atoms will be read.
03910     
03911     // skip lines before number of atoms
03912     for (int i = 0; i < 3; i++) {
03913         fgets(line, lineLen, outputFile); 
03914     }
03915     
03916     iret = fscanf(outputFile,"%d\n", &numAtomsInFile);
03917     if ( iret != 1 ) {
03918         NAMD_die("Error reading QM forces file.");
03919     }
03920     
03921     #ifdef DEBUG_QM
03922     if(numAtomsInFile != msg->numAllAtoms) {
03923         NAMD_die("Error reading QM forces file. Number of atoms in QM output\
03924         does not match the expected.");
03925     }
03926     #endif
03927     
03928     // skip lines before energy
03929     for (int i = 0; i < 3; i++) {
03930         fgets(line, lineLen, outputFile); 
03931     }
03932     
03933     iret = fscanf(outputFile,"%lf\n", &resMsg->energyOrig);
03934     if ( iret != 1 ) {
03935         NAMD_die("Error reading QM forces file.");
03936     }
03937 //     iout << "Energy in step (Hartree): " << resMsg->energyOrig << "\n" <<  endi;
03938     // All energies are given in Eh (Hartree)
03939     // NAMD needs energies in kcal/mol
03940     // The conversion factor is 627.509469
03941     resMsg->energyOrig *= 627.509469;
03942 //     iout << "Energy in step (Kcal/mol): " << resMsg->energyOrig << "\n" <<  endi;
03943     
03944     resMsg->energyCorr = resMsg->energyOrig;
03945     
03946     // skip lines before gradient
03947     for (int i = 0; i < 3; i++) {
03948         fgets(line, lineLen, outputFile) ;
03949     }
03950     
03951     // Break the loop when we find all gradients for QM atoms, 
03952     // and keep the next gradient lines from being read, if there
03953     // are more. Following gradients, if present, will refer to link
03954     // atoms.
03955     atmIndx = 0;
03956     gradCount = 0;
03957     for ( size_t i=0; i<msg->numAllAtoms*3; ++i ) {
03958         
03959         iret = fscanf(outputFile,"%lf\n", &gradient[gradCount]);
03960         if ( iret != 1 ) { NAMD_die("Error reading QM forces file."); }
03961         
03962         if (gradCount == 2){
03963             
03964             // All gradients are given in Eh/a0 (Hartree over Bohr radius)
03965             // NAMD needs forces in kcal/mol/angstrons
03966             // The conversion factor is 627.509469/0.529177 = 1185.82151
03967             if (atmIndx < msg->numQMAtoms ) {
03968                 resForce[atmIndx].force.x = -1*gradient[0]*1185.82151;
03969                 resForce[atmIndx].force.y = -1*gradient[1]*1185.82151;
03970                 resForce[atmIndx].force.z = -1*gradient[2]*1185.82151;
03971             }
03972             
03973             // If we are reading forces applied on Dummy atoms,
03974             // place them on the appropriate QM and MM atoms to conserve energy.
03975             
03976             // This implementation was based on the description in 
03977             // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
03978             if ( msg->numQMAtoms <= atmIndx &&
03979             atmIndx < msg->numAllAtoms ) {
03980                 // The dummy atom points to the QM atom to which it is bound.
03981                 // The QM atom and the MM atom (in a QM-MM bond) point to each other.
03982                 int qmInd = atmP[atmIndx].bountToIndx ;
03983                 int mmInd = atmP[qmInd].bountToIndx ;
03984                 
03985                 Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
03986                 Real mmqmDist = dir.length() ;
03987                 
03988                 Real linkDist = Vector(atmP[atmIndx].position - atmP[qmInd].position).length() ;
03989                 
03990                 Force mmForce(0), qmForce(0), linkForce(gradient[0], gradient[1], gradient[2]);
03991                 linkForce *= -1*1185.82151;
03992                 
03993                 Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
03994                 // Unit vectors
03995                 Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
03996                 Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
03997                 Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
03998                 Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
03999                 
04000                 qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv + 
04001                             (xDelta)*base) )*xuv;
04002                 
04003                 qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv + 
04004                             (yDelta)*base) )*yuv;
04005                 
04006                 qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv + 
04007                             (zDelta)*base) )*zuv;
04008                 
04009                 
04010                 mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
04011                             (xDelta)*base) )*xuv;
04012                 
04013                 mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
04014                             (yDelta)*base) )*yuv;
04015                 
04016                 mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
04017                             (zDelta)*base) )*zuv;
04018                 
04019                 resForce[qmInd].force += qmForce;
04020                 resForce[msg->numQMAtoms + mmInd].force += mmForce;
04021             }
04022             
04023             gradCount = 0;
04024             atmIndx++ ;
04025         } else
04026             gradCount++ ;
04027         
04028     }
04029     
04030     fclose(outputFile);
04031     
04032     // In case charges are not to be read form the QM software output,
04033     // we load the origianl atom charges.
04034     if (msg->qmAtmChrgMode == QMCHRGNONE) {
04035         int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
04036         const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
04037         Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
04038         
04039         atmIndx = 0 ;
04040         for (; atmIndx < msg->numQMAtoms; atmIndx++) {
04041             
04042             // Loops over all QM atoms (in all QM groups) comparing their global indices
04043             for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
04044                 
04045                 if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
04046                     
04047                     atmP[atmIndx].charge = qmAtmChrg[qmIter];
04048                     resForce[atmIndx].charge = qmAtmChrg[qmIter];
04049                     
04050                     break;
04051                 }
04052                 
04053             }
04054             
04055         }
04056     }
04057     else {
04058         // opens redirected output file
04059         outputFile = fopen(tmpRedirectFileName.c_str(),"r");
04060         if ( ! outputFile ) {
04061             iout << iERROR << "Could not find Redirect output file:"
04062             << tmpRedirectFileName << std::endl << endi;
04063             NAMD_err(0); 
04064         }
04065         
04066         DebugM(4,"Opened tmeporary output for charge reading: " << tmpRedirectFileName.c_str() << "\n");
04067         
04068         int lineState = 0;
04069         char result[20] ;
04070         while ( fgets(line, lineLen, outputFile) != NULL){
04071             
04072             // We loop over the lines of the output file untill we find
04073             // the line that initiates the charges lines. Then
04074             // we read all lines and expect to get one charge
04075             // per line, untill we find a line that sums all charges.
04076             
04077             switch (msg->qmAtmChrgMode) {
04078                 
04079                 case QMCHRGMULLIKEN:
04080                 {
04081                 
04082                 if ( strstr(line,"MULLIKEN ATOMIC CHARGES") != NULL ) {
04083                     lineState = 1;
04084                     atmIndx = 0;
04085                     
04086                     // Now we expect the following line to have a series of dashes
04087                     // and the folowing lines to have atom charges (among other info)
04088                     continue;
04089                 }
04090                 
04091                 if ( strstr(line,"Sum of atomic charges") != NULL ) {
04092                     
04093                     // Now that all charges have been read, grabs the total charge
04094                     // just for fun... and checking purposes.
04095                     #ifdef DEBUG_QM
04096                     strncpy(result, line + 31,12) ;
04097                     result[12] = '\0';
04098                     
04099                     DebugM(4,"Total charge of QM region calculated by ORCA is: " 
04100                     << atof(result) << std::endl )
04101                     #endif
04102                     
04103                     // Nothing more to read from the output file
04104                     break;
04105                 }
04106                 
04107                 // Line state 1 means we have to skip ONLY one line.
04108                 if (lineState == 1) {
04109                     lineState = 2;
04110                     continue;
04111                 }
04112                 
04113                 // Line state 2 means we have to read the line and grab the info.
04114                 if (lineState == 2) {
04115                     
04116                     strncpy(result, line+8,12) ;
04117                     result[12] = '\0';
04118                     
04119                     Real localCharge = atof(result);
04120                     
04121                     // If we are reading charges from QM atoms, store them.
04122                     if (atmIndx < msg->numQMAtoms ) {
04123                         atmP[atmIndx].charge = localCharge;
04124                         resForce[atmIndx].charge = localCharge;
04125                     }
04126                     
04127                     // If we are reading charges from Dummy atoms,
04128                     // place the on the appropriate QM atom.
04129                      if ( msg->numQMAtoms <= atmIndx &&
04130                         atmIndx < msg->numAllAtoms ) {
04131                         int qmInd = atmP[atmIndx].bountToIndx ;
04132                         atmP[qmInd].charge += localCharge;
04133                         resForce[qmInd].charge += localCharge;
04134                     }
04135                     
04136                     atmIndx++ ;
04137                     
04138                     // If we found all charges for QM atoms, change the lineState
04139                     // untill we reach the "total charge" line.
04140                     if (atmIndx == msg->numAllAtoms ) {
04141                         lineState = 0;
04142                     }
04143                     
04144                     continue;
04145                 }
04146                 
04147                 } break ;
04148                 
04149                 case QMCHRGCHELPG :
04150                 {
04151                 
04152                 if ( strstr(line,"CHELPG Charges") != NULL ) {
04153                     lineState = 1;
04154                     atmIndx = 0;
04155                     
04156                     // Now we expect the following line to have a series of dashes
04157                     // and the folowing lines to have atom charges (among other info)
04158                     continue;
04159                 }
04160                 
04161                 if ( strstr(line,"Total charge") != NULL ) {
04162                     
04163                     // Now that all charges have been read, grabs the total charge
04164                     // just for fun... and checking purposes.
04165                     #ifdef DEBUG_QM
04166                     strncpy(result, line + 14,13) ;
04167                     result[13] = '\0';
04168                     
04169                     DebugM(4,"Total charge of QM region calculated by ORCA is: " 
04170                     << atof(result) << std::endl )
04171                     #endif
04172                     
04173                     // Nothing more to read from the output file
04174                     break;
04175                 }
04176                 
04177                 // Line state 1 means we have to skip ONLY one line.
04178                 if (lineState == 1) {
04179                     lineState = 2;
04180                     continue;
04181                 }
04182                 
04183                 // Line state 2 means we have to read the line and grab the info.
04184                 if (lineState == 2) {
04185                     
04186                     strncpy(result, line+12,15) ;
04187                     result[15] = '\0';
04188                     
04189                     Real localCharge = atof(result);
04190                     
04191                     // If we are reading charges from QM atoms, store them.
04192                     if (atmIndx < msg->numQMAtoms ) {
04193                         atmP[atmIndx].charge = localCharge;
04194                         resForce[atmIndx].charge = localCharge;
04195                     }
04196                     
04197                     // If we are reading charges from Dummy atoms,
04198                     // place the on the appropriate QM atom.
04199                      if ( msg->numQMAtoms <= atmIndx &&
04200                         atmIndx < msg->numAllAtoms ) {
04201                         int qmInd = atmP[atmIndx].bountToIndx ;
04202                         atmP[qmInd].charge += localCharge;
04203                         resForce[qmInd].charge += localCharge;
04204                     }
04205                     
04206                     atmIndx++ ;
04207                     
04208                     // If we found all charges for QM atoms, we ignore the following line
04209                     // untill we reach the "total charge" line.
04210                     if (atmIndx == msg->numAllAtoms ) {
04211                         lineState = 1;
04212                     }
04213                     
04214                     continue;
04215                 }
04216                 
04217                 } break;
04218             }
04219         }
04220         
04221         fclose(outputFile);
04222     }
04223     
04224     // remove force file
04225 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
04226 //     iret = remove(outputFileName);
04227 //     if ( iret ) { NAMD_err(0); }
04228     
04229     
04230     DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
04231     
04232     atmP = msg->data ;
04233     pcP = msg->data + msg->numAllAtoms ;
04234     
04235     // The initial point charge index for the force message is the number of QM
04236     // atoms, since the dummy atoms have no representation in NAMD
04237     int pcIndx = msg->numQMAtoms;
04238     
04239     // If there are no point charges, orca will not create a "pcgrad" file.
04240     if (msg->numAllPntChrgs > 0) {
04241     
04242         // opens redirected output file
04243         outputFile = fopen(pcGradFileName.c_str(),"r");
04244         if ( ! outputFile ) {
04245             iout << iERROR << "Could not find PC gradient output file:"
04246             << pcGradFileName << std::endl << endi;
04247             NAMD_err(0); 
04248         }
04249         
04250         DebugM(4,"Opened pc output for gradient reading: " << pcGradFileName.c_str() << "\n");
04251         
04252         int pcgradNumPC = 0, readPC = 0;
04253         iret = fscanf(outputFile,"%d\n", &pcgradNumPC );
04254         if ( iret != 1 ) {
04255             NAMD_die("Error reading PC forces file.");
04256         }
04257         DebugM(4,"Found in pc gradient output: " << pcgradNumPC << " point charge grads.\n");
04258         
04259         // We loop over point charges, reading the total electrostatic force 
04260         // applied on them by the QM region.
04261         // We redistribute the forces applied over virtual point
04262         // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
04263         for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
04264             
04265             Force totalForce(0);
04266             
04267             // No force was applied to the QM region due to this charge, since it
04268             // was skipped when writing down point charges to the QM software, so it
04269             // does not receive any back from the QM region. It must be an MM1 atom
04270             // from a QM-MM bond.
04271             if (pcP[i].type == QMPCTYPE_IGNORE)
04272                 continue;
04273             
04274             fgets(line, lineLen, outputFile);
04275             
04276             iret = sscanf(line,"%lf %lf %lf\n", &totalForce[0], &totalForce[1], &totalForce[2] );
04277             if ( iret != 3 ) {
04278                 NAMD_die("Error reading PC forces file.");
04279             }
04280             // All gradients are given in Eh/a0 (Hartree over Bohr radius)
04281             // NAMD needs forces in kcal/mol/angstrons
04282             // The conversion factor is 627.509469/0.529177 = 1185.82151
04283             totalForce *= -1185.82151;
04284             
04285             if (pcP[i].type == QMPCTYPE_CLASSICAL) {
04286                 // If we already ignored MM1 charges, we take all other 
04287                 // non-virtual charges and apply forces directly to them.
04288                 resForce[pcIndx].force += totalForce;
04289             }
04290             else {
04291                 // If we are handling virtual PC, we distribute the force over
04292                 // MM1 and MM2.
04293                 
04294                 Force mm1Force(0), mm2Force(0);
04295                 
04296                 // Virtual PC are bound to MM2.
04297                 int mm2Indx = pcP[i].bountToIndx;
04298                 // MM2 charges are bound to MM1.
04299                 int mm1Indx = pcP[mm2Indx].bountToIndx;
04300                 
04301                 Real Cq = pcP[i].dist;
04302                 
04303                 mm1Force = (1-Cq)*totalForce ;
04304                 mm2Force = Cq*totalForce ;
04305                 
04306                 resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
04307                 resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
04308             }
04309             
04310         }
04311         
04312         fclose(outputFile);
04313     }
04314     
04315     delete [] line;
04316     
04317     // Adjusts forces from PME, canceling contributions from the QM and 
04318     // direct Coulomb forces calculated here.
04319     if (msg->PMEOn) {
04320         
04321         DebugM(1,"Correcting forces and energy for PME.\n");
04322         
04323         ewaldcof = msg->PMEEwaldCoefficient;
04324         BigReal TwoBySqrtPi = 1.12837916709551;
04325         pi_ewaldcof = TwoBySqrtPi * ewaldcof;
04326         
04327         for (size_t i=0; i < msg->numQMAtoms; i++) {
04328             
04329             BigReal p_i_charge = atmP[i].charge ;
04330             Position pos_i = atmP[i].position ;
04331             
04332             for (size_t j=i+1; j < msg->numQMAtoms; j++) {
04333                 
04334                 BigReal p_j_charge = atmP[j].charge ;
04335                 
04336                 Position pos_j = atmP[j].position ;
04337                 
04338                 BigReal r = Vector(pos_i - pos_j).length() ;
04339                 
04340                 BigReal tmp_a = r * ewaldcof;
04341                 BigReal tmp_b = erfc(tmp_a);
04342                 BigReal corr_energy = tmp_b;
04343                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
04344                 
04345 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
04346                 BigReal recip_energy = (1-tmp_b)/r;
04347                 
04348 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
04349                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
04350                 
04351                 const BigReal kq_i = p_i_charge * constants;
04352                 
04353                 // Final force and energy correction for this pair of atoms.
04354                 BigReal energy = kq_i * p_j_charge * recip_energy ;
04355                 
04356                 Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
04357                 
04358                 // The force is *subtracted* from the total force acting on
04359                 // both atoms. The sign on fixForce corrects the orientation
04360                 // of the subtracted force.
04361 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
04362 //                     << std::endl);
04363 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
04364 //                     << std::endl);
04365 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
04366                 resForce[i].force -= fixForce ;
04367                 resForce[j].force -= -1*fixForce;
04368                 
04369                 // The energy is *subtracted* from the total energy calculated here.
04370                 resMsg->energyCorr -= energy;
04371             }
04372         }
04373         
04374         pcIndx = msg->numQMAtoms;
04375         // We only loop over point charges from real atoms, ignoring the ones 
04376         // created to handle QM-MM bonds.
04377         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
04378             
04379             BigReal p_i_charge = pcPpme[i].charge;
04380             Position pos_i = pcPpme[i].position ;
04381             
04382             const BigReal kq_i = p_i_charge * constants;
04383             
04384             Force fixForce = 0;
04385             
04386             for (size_t j=0; j<msg->numQMAtoms; ++j ) {
04387                 
04388                 BigReal p_j_charge = atmP[j].charge ;
04389                 
04390                 Position pos_j = atmP[j].position ;
04391                 
04392                 BigReal r = Vector(pos_i - pos_j).length() ;
04393                 
04394                 BigReal tmp_a = r * ewaldcof;
04395                 BigReal tmp_b = erfc(tmp_a);
04396                 BigReal corr_energy = tmp_b;
04397                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
04398                 
04399 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
04400                 BigReal recip_energy = (1-tmp_b)/r;
04401                 
04402 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
04403                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
04404                 
04405                 // Final force and energy correction for this pair of atoms.
04406                 BigReal energy = kq_i * p_j_charge * recip_energy ;
04407                 
04408                 fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
04409                 
04410                 // The energy is *subtracted* from the total energy calculated here.
04411                 resMsg->energyCorr -= energy;
04412                 
04413             }
04414             
04415             // The force is *subtracted* from the total force acting on
04416                 // the point charge.
04417 //                 DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
04418 //                     << std::endl);
04419 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
04420             resForce[pcIndx].force -= kq_i*fixForce ;
04421         }
04422         
04423     }
04424     
04425     DebugM(1,"Determining virial...\n");
04426     
04427     // Calculates the virial contribution form the forces on QM atoms and 
04428     // point charges calculated here.
04429     for (size_t i=0; i < msg->numQMAtoms; i++) {
04430         // virial used by NAMD is -'ve of normal convention, so reverse it!
04431         // virial[i][j] in file should be sum of -1 * f_i * r_j
04432         for ( int k=0; k<3; ++k )
04433             for ( int l=0; l<3; ++l )
04434                 resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
04435     }
04436     
04437     pcIndx = msg->numQMAtoms; // Index in the real PC force array.
04438     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
04439         // virial used by NAMD is -'ve of normal convention, so reverse it!
04440         // virial[i][j] in file should be sum of -1 * f_i * r_j
04441         for ( int k=0; k<3; ++k )
04442             for ( int l=0; l<3; ++l )
04443                 resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
04444     }
04445     
04446     DebugM(1,"End of QM processing. Sending result message.\n");
04447     
04448     // Send message to rank zero with results.
04449     QMProxy[0].recvQMRes(resMsg);
04450     
04451     if (msg->switching && msg->PMEOn)
04452         delete [] pcPpme;
04453     
04454     delete msg;
04455     return ;
04456 }

void ComputeQMMgr::calcUSR ( QMGrpCalcMsg  ) 

Definition at line 4458 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMForce::charge, QMAtomData::charge, QMGrpCalcMsg::constants, QMGrpCalcMsg::data, DebugM, endi(), QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, QMGrpCalcMsg::execPath, QMGrpResMsg::force, QMGrpCalcMsg::grpIndx, QMGrpResMsg::grpIndx, QMAtomData::id, iERROR(), if(), iout, j, Vector::length(), NAMD_die(), NAMD_err(), QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMGrpResMsg::numForces, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, QMGrpCalcMsg::peIter, QMGrpCalcMsg::PMEOn, QMAtomData::position, QMGrpCalcMsg::prepProc, QMGrpCalcMsg::prepProcOn, QMPCTYPE_IGNORE, QMGrpCalcMsg::secProc, QMGrpCalcMsg::secProcOn, QMGrpCalcMsg::switching, QMGrpCalcMsg::timestep, QMAtomData::type, QMGrpResMsg::virial, Vector::x, x, Vector::y, y, Vector::z, and z.

04458                                             {
04459     
04460     FILE *inputFile,*outputFile;
04461     int iret;
04462     
04463     std::string qmCommand, inputFileName, outputFileName;
04464     
04465     // For coulomb calculation
04466     BigReal constants = msg->constants;
04467     
04468     DebugM(4,"Running USER DEFINED SOFTWARE on PE " << CkMyPe() << std::endl);
04469     
04470     QMAtomData *pcP = msg->data + msg->numAllAtoms ;
04471     
04472     // We create a pointer for PME correction, which may point to
04473     // a copy of the original point charge array, with unchanged charges,
04474     // or points to the original array in case no switching or charge
04475     // scheme is being used.
04476     QMAtomData *pcPpme = NULL;
04477     if (msg->switching) {
04478         
04479         if (msg->PMEOn)
04480             pcPpme = new QMAtomData[msg->numRealPntChrgs];
04481         
04482         pntChrgSwitching(msg, pcPpme) ;
04483     } else {
04484         pcPpme = pcP;
04485     }
04486     
04487     // For each QM group, create a subdirectory where files will be palced.
04488     std::string baseDir(msg->baseDir);
04489     baseDir.append("/") ;
04490     std::ostringstream itosConv ;
04491     itosConv << msg->peIter ;
04492     baseDir += itosConv.str() ;
04493     
04494     struct stat info;
04495     
04496     if (stat(msg->baseDir, &info) != 0 ) {
04497         CkPrintf( "Node %d cannot access directory %s\n",
04498                   CkMyPe(), baseDir.c_str() );
04499         NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
04500     }
04501     else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
04502         DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
04503         int retVal = mkdir(baseDir.c_str(), S_IRWXU);
04504     }
04505     
04506     #ifdef DEBUG_QM
04507     Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
04508     #endif
04509     
04510     inputFileName.clear();
04511     inputFileName.append(baseDir.c_str()) ;
04512     inputFileName.append("/qmmm_") ;
04513     inputFileName += itosConv.str() ;
04514     inputFileName.append(".input") ;
04515     
04516     // Opens file for coordinate and parameter input
04517     inputFile = fopen(inputFileName.c_str(),"w");
04518     if ( ! inputFile ) {
04519         iout << iERROR << "Could not open input file for writing: " 
04520         << inputFileName << "\n" << endi ;
04521         NAMD_err("Error writing QM input file.");
04522     }
04523     
04524     // Builds the command that will be executed
04525     qmCommand.clear();
04526     qmCommand.append("cd ");
04527     qmCommand.append(baseDir);
04528     qmCommand.append(" ; ");
04529     qmCommand.append(msg->execPath) ;
04530     qmCommand.append(" ") ;
04531     qmCommand.append(inputFileName) ;
04532     
04533     // Builds the file name where orca will place the gradient
04534     // This will be relative to the input file
04535     outputFileName = inputFileName ;
04536     outputFileName.append(".result") ;
04537     
04538     int numPntChrgs = 0;
04539     for (int i=0; i<msg->numAllPntChrgs; i++ ) {
04540         if (pcP[i].type != QMPCTYPE_IGNORE)
04541             numPntChrgs++;
04542     }
04543     
04544     iret = fprintf(inputFile,"%d %d\n",msg->numAllAtoms, numPntChrgs);
04545     if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
04546     
04547     DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " << 
04548         inputFileName.c_str() << " and " << msg->numAllPntChrgs << 
04549         " point charges." << std::endl);
04550     
04551     // write QM and dummy atom coordinates to input file.
04552     QMAtomData *atmP = msg->data ;
04553     for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
04554         
04555         double x = atmP->position.x;
04556         double y = atmP->position.y;
04557         double z = atmP->position.z;
04558         
04559         iret = fprintf(inputFile,"%f %f %f %s\n",
04560                        x,y,z,atmP->element);
04561         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
04562         
04563     }
04564     
04565     int numWritenPCs = 0;
04566     // Write point charges to file.
04567     pcP = msg->data + msg->numAllAtoms ;
04568     for ( size_t j=0; j < msg->numAllPntChrgs; j++, ++pcP) {
04569         
04570         if (pcP->type == QMPCTYPE_IGNORE)
04571                 continue;
04572         
04573         double charge = pcP->charge;
04574         
04575         double x = pcP->position.x;
04576         double y = pcP->position.y;
04577         double z = pcP->position.z;
04578         
04579         iret = fprintf(inputFile,"%f %f %f %f\n",
04580                        x,y,z,charge);
04581         if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
04582         
04583         numWritenPCs++;
04584     }
04585     
04586     DebugM(4,"Closing input file\n");
04587     fclose(inputFile);
04588     
04589     if (msg->prepProcOn) {
04590         
04591         std::string prepProc(msg->prepProc) ;
04592         prepProc.append(" ") ;
04593         prepProc.append(inputFileName) ;
04594         iret = system(prepProc.c_str());
04595         if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
04596         if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
04597     }
04598     
04599         // runs QM command
04600     DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
04601     iret = system(qmCommand.c_str());
04602     
04603     if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
04604     if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
04605 
04606     if (msg->secProcOn) {
04607         
04608         std::string secProc(msg->secProc) ;
04609         secProc.append(" ") ;
04610         secProc.append(inputFileName) ;
04611         itosConv.str("");
04612         itosConv.clear() ;
04613         itosConv << msg->timestep ;
04614         secProc.append(" ") ;
04615         secProc += itosConv.str() ;
04616         
04617         iret = system(secProc.c_str());
04618         if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
04619         if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
04620     }
04621 
04622     // remove coordinate file
04623 //     iret = remove(inputFileName);
04624 //     if ( iret ) { NAMD_err(0); }
04625 
04626     // remove coordinate file
04627 //     iret = remove(pntChrgFileName);
04628 //     if ( iret ) { NAMD_err(0); }
04629     
04630     // opens output file
04631     DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
04632     outputFile = fopen(outputFileName.c_str(),"r");
04633     if ( ! outputFile ) {
04634         iout << iERROR << "Could not find QM output file!\n" << endi;
04635         NAMD_err(0); 
04636     }
04637 
04638     // Resets the pointers.
04639     atmP = msg->data ;
04640     pcP = msg->data + msg->numAllAtoms ;
04641     
04642     // Allocates the return message.
04643     QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
04644     resMsg->grpIndx = msg->grpIndx;
04645     resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
04646     resMsg->energyOrig = 0;
04647     resMsg->energyCorr = 0;
04648     for ( int k=0; k<3; ++k )
04649         for ( int l=0; l<3; ++l )
04650             resMsg->virial[k][l] = 0;
04651     QMForce *resForce = resMsg->force ;
04652     
04653     // We split the data into two pointers so we can "skip" the dummy atoms
04654     // which have no representation in NAMD.
04655     for (int i=0; i<resMsg->numForces; i++) {
04656         resForce[i].force = 0;
04657         resForce[i].charge = 0 ;
04658         if (i < msg->numQMAtoms) {
04659             // We use the replace field to indicate QM atoms,
04660             // which will have charge information.
04661             resForce[i].replace = 1 ;
04662             resForce[i].id = atmP->id;
04663             atmP++;
04664         }
04665         else {
04666             // We use the replace field to indicate QM atoms,
04667             // which will have charge information.
04668             resForce[i].replace = 0 ;
04669             resForce[i].id = pcP->id;
04670             pcP++;
04671         }
04672     }
04673     
04674     // Resets the pointers.
04675     atmP = msg->data ;
04676     pcP = msg->data + msg->numAllAtoms ;
04677     
04678     // Reads the data form the output file created by the QM software.
04679     // Gradients over the QM atoms, and Charges for QM atoms will be read.
04680     
04681     // Number of point charges for which we will receive forces.
04682     int usrPCnum = 0;
04683     const size_t lineLen = 256;
04684     char *line = new char[lineLen];
04685     
04686     fgets(line, lineLen, outputFile);
04687     
04688 //     iret = fscanf(outputFile,"%lf %d\n", &resMsg->energyOrig, &usrPCnum);
04689     iret = sscanf(line,"%lf %i\n", &resMsg->energyOrig, &usrPCnum);
04690     if ( iret < 1 ) {
04691         NAMD_die("Error reading energy from QM results file.");
04692     }
04693     
04694     resMsg->energyCorr = resMsg->energyOrig;
04695     
04696     if (iret == 2 && numWritenPCs != usrPCnum) {
04697         iout << iERROR << "Number of point charges does not match what was provided!\n" << endi ;
04698         NAMD_die("Error reading QM results file.");
04699     }
04700     
04701     size_t atmIndx;
04702     double localForce[3];
04703     double localCharge;
04704     for (atmIndx = 0; atmIndx < msg->numAllAtoms; atmIndx++) {
04705         
04706         iret = fscanf(outputFile,"%lf %lf %lf %lf\n", 
04707                       localForce+0,
04708                       localForce+1,
04709                       localForce+2,
04710                       &localCharge);
04711         if ( iret != 4 ) {
04712             NAMD_die("Error reading forces and charge from QM results file.");
04713         }
04714         
04715         // If we are reading charges and forces on QM atoms, store
04716         // them directly.
04717         if (atmIndx < msg->numQMAtoms ) {
04718             
04719             resForce[atmIndx].force.x = localForce[0];
04720             resForce[atmIndx].force.y = localForce[1];
04721             resForce[atmIndx].force.z = localForce[2];
04722             
04723             atmP[atmIndx].charge = localCharge;
04724             resForce[atmIndx].charge = localCharge;
04725         }
04726         
04727         // If we are reading charges and forces applied on Dummy atoms,
04728         // place them on the appropriate QM and MM atoms to conserve energy.
04729         
04730         // This implementation was based on the description in 
04731         // DOI: 10.1002/jcc.20857  :  Equations 30 to 32
04732         if ( msg->numQMAtoms <= atmIndx &&
04733         atmIndx < msg->numAllAtoms ) {
04734             
04735             // We keep the calculated charge in the dummy atom
04736             // so that we can calculate electrostatic forces later on.
04737             atmP[atmIndx].charge = localCharge;
04738             
04739             // If we are reading charges from Dummy atoms,
04740             // place them on the appropriate QM atom.
04741             // The dummy atom points to the QM atom to which it is bound.
04742             int qmInd = atmP[atmIndx].bountToIndx ;
04743             resForce[qmInd].charge += localCharge;
04744             
04745             // The dummy atom points to the QM atom to which it is bound.
04746             // The QM atom and the MM atom (in a QM-MM bond) point to each other.
04747             int mmInd = atmP[qmInd].bountToIndx ;
04748             
04749             Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
04750             Real mmqmDist = dir.length() ;
04751             
04752             Real linkDist = Vector(atmP[atmIndx].position - 
04753                             atmP[qmInd].position).length() ;
04754             
04755             Force mmForce(0), qmForce(0), 
04756                 linkForce(localForce[0], localForce[1], localForce[2]);
04757             
04758             Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
04759             // Unit vectors
04760             Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
04761             Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
04762             Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
04763             Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
04764             
04765             qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv + 
04766                         (xDelta)*base) )*xuv;
04767             
04768             qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv + 
04769                         (yDelta)*base) )*yuv;
04770             
04771             qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv + 
04772                         (zDelta)*base) )*zuv;
04773             
04774             
04775             mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
04776                         (xDelta)*base) )*xuv;
04777             
04778             mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
04779                         (yDelta)*base) )*yuv;
04780             
04781             mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
04782                         (zDelta)*base) )*zuv;
04783             
04784             resForce[qmInd].force += qmForce;
04785             resForce[msg->numQMAtoms + mmInd].force += mmForce;
04786         }
04787     }
04788     
04789     // The initial point charge index for the force message is the number of QM
04790     // atoms, since the dummy atoms have no representation in NAMD
04791     int pcIndx = msg->numQMAtoms;
04792     
04793     if (usrPCnum > 0) {
04794         // We loop over point charges, reading the total electrostatic force 
04795         // applied on them by the QM region.
04796         // We redistribute the forces applied over virtual point
04797         // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
04798         for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
04799             
04800             Force totalForce(0);
04801             
04802             // No force was applied to the QM region due to this charge, since it
04803             // was skipped when writing down point charges to the QM software, so it
04804             // does not receive any back from the QM region. It must be an MM1 atom
04805             // from a QM-MM bond.
04806             if (pcP[i].type == QMPCTYPE_IGNORE)
04807                 continue;
04808             
04809             iret = fscanf(outputFile,"%lf %lf %lf\n", 
04810                            &totalForce[0], &totalForce[1], &totalForce[2]);
04811             if ( iret != 3 ) {
04812                 NAMD_die("Error reading PC forces from QM results file.");
04813             }
04814             
04815             if (pcP[i].type == QMPCTYPE_CLASSICAL) {
04816                 // If we already ignored MM1 charges, we take all other 
04817                 // non-virtual charges and apply forces directly to them.
04818                 resForce[pcIndx].force += totalForce;
04819             }
04820             else {
04821                 // If we are handling virtual PC, we distribute the force over
04822                 // MM1 and MM2.
04823                 
04824                 Force mm1Force(0), mm2Force(0);
04825                 
04826                 // Virtual PC are bound to MM2.
04827                 int mm2Indx = pcP[i].bountToIndx;
04828                 // MM2 charges are bound to MM1.
04829                 int mm1Indx = pcP[mm2Indx].bountToIndx;
04830                 
04831                 Real Cq = pcP[i].dist;
04832                 
04833                 mm1Force = (1-Cq)*totalForce ;
04834                 mm2Force = Cq*totalForce ;
04835                 
04836                 resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
04837                 resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
04838             }
04839             
04840             
04841         }
04842     }
04843     
04844     fclose(outputFile);
04845     delete [] line;
04846     
04847     // In case charges are not to be read form the QM software output,
04848     // we load the origianl atom charges.
04849     if (msg->qmAtmChrgMode == QMCHRGNONE) {
04850         int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
04851         const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
04852         Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
04853         
04854         atmIndx = 0 ;
04855         for (; atmIndx < msg->numQMAtoms; atmIndx++) {
04856             
04857             // Loops over all QM atoms (in all QM groups) comparing their global indices
04858             for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
04859                 
04860                 if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
04861                     
04862                     atmP[atmIndx].charge = qmAtmChrg[qmIter];
04863                     resForce[atmIndx].charge = qmAtmChrg[qmIter];
04864                     
04865                     break;
04866                 }
04867                 
04868             }
04869             
04870         }
04871         for (size_t j=msg->numQMAtoms; j<msg->numAllAtoms; ++j ) {
04872             atmP[j].charge = 0;
04873         }
04874     }
04875     
04876     // remove force file
04877 //     DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
04878 //     iret = remove(outputFileName);
04879 //     if ( iret ) { NAMD_err(0); }
04880     
04881     if (usrPCnum == 0) {
04882         DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
04883         
04884         atmP = msg->data ;
04885         pcP = msg->data + msg->numAllAtoms ;
04886         
04887         // We only loop over point charges from real atoms, ignoring the ones 
04888         // created to handle QM-MM bonds.
04889         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
04890             
04891             // No force was applied to the QM region due to this charge, so it 
04892             // does not receive any back from the QM region. It must be an MM1 atom
04893             // from a QM-MM bond.
04894             if (pcP[i].type == QMPCTYPE_IGNORE)
04895                 continue;
04896             
04897             Force totalForce(0);
04898             
04899             BigReal pntCharge = pcP[i].charge;
04900             
04901             Position posMM = pcP[i].position ;
04902             
04903             for (size_t j=0; j<msg->numAllAtoms; ++j ) {
04904                 
04905                 BigReal qmCharge = atmP[j].charge ;
04906                 
04907                 BigReal force = pntCharge*qmCharge*constants ;
04908                 
04909                 Position rVec = posMM - atmP[j].position ;
04910                 
04911                 force /= rVec.length2();
04912                 
04913                 // We accumulate the total force felt by a point charge
04914                 // due to the QM charge distribution. This total force
04915                 // will be applied to the point charge if it is a real one,
04916                 // or will be distirbuted over MM1 and MM2 point charges, it 
04917                 // this is a virtual point charge.
04918                 totalForce += force*rVec.unit();
04919             }
04920             
04921             if (pcP[i].type == QMPCTYPE_CLASSICAL) {
04922                 // If we already ignored MM1 charges, we take all other 
04923                 // non-virtual charges and apply forces directly to them.
04924                 resForce[pcIndx].force += totalForce;
04925             }
04926             else {
04927                 // If we are handling virtual PC, we distribute the force over
04928                 // MM1 and MM2.
04929                 
04930                 Force mm1Force(0), mm2Force(0);
04931                 
04932                 // Virtual PC are bound to MM2.
04933                 int mm2Indx = pcP[i].bountToIndx;
04934                 // MM2 charges are bound to MM1.
04935                 int mm1Indx = pcP[mm2Indx].bountToIndx;
04936                 
04937                 Real Cq = pcP[i].dist;
04938                 
04939                 mm1Force = (1-Cq)*totalForce ;
04940                 mm2Force = Cq*totalForce ;
04941                 
04942                 resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
04943                 resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
04944             }
04945             
04946         }
04947     }
04948     
04949     // Adjusts forces from PME, canceling contributions from the QM and 
04950     // direct Coulomb forces calculated here.
04951     if (msg->PMEOn) {
04952         
04953         DebugM(1,"Correcting forces and energy for PME.\n");
04954         
04955         ewaldcof = msg->PMEEwaldCoefficient;
04956         BigReal TwoBySqrtPi = 1.12837916709551;
04957         pi_ewaldcof = TwoBySqrtPi * ewaldcof;
04958         
04959         for (size_t i=0; i < msg->numQMAtoms; i++) {
04960             
04961             BigReal p_i_charge = atmP[i].charge ;
04962             Position pos_i = atmP[i].position ;
04963             
04964             for (size_t j=i+1; j < msg->numQMAtoms; j++) {
04965                 
04966                 BigReal p_j_charge = atmP[j].charge ;
04967                 
04968                 Position pos_j = atmP[j].position ;
04969                 
04970                 BigReal r = Vector(pos_i - pos_j).length() ;
04971                 
04972                 BigReal tmp_a = r * ewaldcof;
04973                 BigReal tmp_b = erfc(tmp_a);
04974                 BigReal corr_energy = tmp_b;
04975                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
04976                 
04977 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
04978                 BigReal recip_energy = (1-tmp_b)/r;
04979                 
04980 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
04981                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
04982                 
04983                 const BigReal kq_i = p_i_charge * constants;
04984                 
04985                 // Final force and energy correction for this pair of atoms.
04986                 BigReal energy = kq_i * p_j_charge * recip_energy ;
04987                 
04988                 Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
04989                 
04990                 // The force is *subtracted* from the total force acting on
04991                 // both atoms. The sign on fixForce corrects the orientation
04992                 // of the subtracted force.
04993 //                 DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
04994 //                     << std::endl);
04995 //                 DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
04996 //                     << std::endl);
04997 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
04998                 resForce[i].force -= fixForce ;
04999                 resForce[j].force -= -1*fixForce;
05000                 
05001                 // The energy is *subtracted* from the total energy calculated here.
05002                 resMsg->energyCorr -= energy;
05003             }
05004         }
05005         
05006         pcIndx = msg->numQMAtoms;
05007         // We only loop over point charges from real atoms, ignoring the ones 
05008         // created to handle QM-MM bonds.
05009         for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
05010             
05011             BigReal p_i_charge = pcPpme[i].charge;
05012             Position pos_i = pcPpme[i].position ;
05013             
05014             const BigReal kq_i = p_i_charge * constants;
05015             
05016             Force fixForce = 0;
05017             
05018             for (size_t j=0; j<msg->numQMAtoms; ++j ) {
05019                 
05020                 BigReal p_j_charge = atmP[j].charge ;
05021                 
05022                 Position pos_j = atmP[j].position ;
05023                 
05024                 BigReal r = Vector(pos_i - pos_j).length() ;
05025                 
05026                 BigReal tmp_a = r * ewaldcof;
05027                 BigReal tmp_b = erfc(tmp_a);
05028                 BigReal corr_energy = tmp_b;
05029                 BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
05030                 
05031 //                 BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
05032                 BigReal recip_energy = (1-tmp_b)/r;
05033                 
05034 //                 BigReal recip_gradient = -(1-corr_gradient)/(r*2);
05035                 BigReal recip_gradient = -(1-corr_gradient)/(r*r);
05036                 
05037                 // Final force and energy correction for this pair of atoms.
05038                 BigReal energy = kq_i * p_j_charge * recip_energy ;
05039                 
05040                 fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
05041                 
05042                 // The energy is *subtracted* from the total energy calculated here.
05043                 resMsg->energyCorr -= energy;
05044                 
05045             }
05046             
05047             // The force is *subtracted* from the total force acting on
05048                 // the point charge.
05049 //                 DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
05050 //                     << std::endl);
05051 //                 DebugM(4,"Force correction: " << fixForce << std::endl);
05052             resForce[pcIndx].force -= kq_i*fixForce ;
05053         }
05054         
05055     }
05056     
05057     DebugM(1,"Determining virial...\n");
05058     
05059     // Calculates the virial contribution form the forces on QM atoms and 
05060     // point charges calculated here.
05061     for (size_t i=0; i < msg->numQMAtoms; i++) {
05062         // virial used by NAMD is -'ve of normal convention, so reverse it!
05063         // virial[i][j] in file should be sum of -1 * f_i * r_j
05064         for ( int k=0; k<3; ++k )
05065             for ( int l=0; l<3; ++l )
05066                 resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
05067     }
05068     
05069     pcIndx = msg->numQMAtoms; // Index in the real PC force array.
05070     for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
05071         // virial used by NAMD is -'ve of normal convention, so reverse it!
05072         // virial[i][j] in file should be sum of -1 * f_i * r_j
05073         for ( int k=0; k<3; ++k )
05074             for ( int l=0; l<3; ++l )
05075                 resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
05076     }
05077     
05078     DebugM(1,"End of QM processing. Sending result message.\n");
05079     
05080     // Send message to rank zero with results.
05081     QMProxy[0].recvQMRes(resMsg);
05082     
05083     if (msg->switching && msg->PMEOn)
05084         delete [] pcPpme;
05085     
05086     delete msg;
05087     return ;
05088 }

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, DebugM, endi(), QMForceMsg::energy, QMForceMsg::force, Molecule::get_numQMAtoms(), QMForce::id, ComputeQMAtom::id, iINFO(), SortedArray< Elem >::insert(), iout, j, QMForceMsg::lssDat, QMPntChrgMsg::numAtoms, QMCoordMsg::numAtoms, QMForceMsg::numForces, QMForceMsg::numLssDat, SimParameters::PMEOn, QMForceMsg::PMEOn, SimParameters::qmOutFreq, SimParameters::qmPosOutFreq, ResizeArray< Elem >::size(), SortedArray< Elem >::sort(), QMPntChrgMsg::sourceNode, QMCoordMsg::sourceNode, QMForceMsg::virial, write_dcdstep(), x, y, and 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         DebugM(1,"Sending forces and charges to source " << j << std::endl);
02541         
02542         QMCoordMsg *qmmsg = 0;
02543         QMPntChrgMsg *pcmsg = 0 ;
02544         
02545         int totForces = 0;
02546         int sourceNode = -1;
02547         
02548         if (pntChrgCoordMsgs == NULL) {
02549             
02550             qmmsg = qmCoordMsgs[j];
02551             qmCoordMsgs[j] = 0;
02552             
02553             totForces = qmmsg->numAtoms ;
02554             
02555             sourceNode = qmmsg->sourceNode;
02556         }
02557         else {
02558             pcmsg = pntChrgCoordMsgs[j];
02559             pntChrgCoordMsgs[j] = 0;
02560             
02561             sourceNode = pcmsg->sourceNode;
02562             
02563             // Since we receive two different messages from nodes, there is no 
02564             // guarantee the two sets of messages will come in the same order.
02565             // Therefore, we match the messages by comaring their sourceNodes.
02566             for (int aux=0; aux<numSources; aux++) {
02567                 
02568                 if (qmCoordMsgs[aux] == 0)
02569                     continue;
02570                 
02571                 qmmsg = qmCoordMsgs[aux];
02572                 
02573                 if (qmmsg->sourceNode == sourceNode) {
02574                     qmCoordMsgs[aux] = 0;
02575                     break;
02576                 }
02577             }
02578             
02579             DebugM(1,"Building force mesage for rank " 
02580             << pcmsg->sourceNode << std::endl);
02581             
02582             totForces = qmmsg->numAtoms + pcmsg->numAtoms;
02583         }
02584         
02585         QMForceMsg *fmsg = new (totForces, subsArray.size(), 0) QMForceMsg;
02586         fmsg->PMEOn = simParams->PMEOn;
02587         fmsg->numForces = totForces;
02588         fmsg->numLssDat = subsArray.size();
02589         
02590         DebugM(1,"Loading QM forces.\n");
02591         
02592         // This iterator is used in BOTH loops!
02593         int forceIter = 0;
02594         
02595         for ( int i=0; i < qmmsg->numAtoms; ++i ) {
02596             fmsg->force[forceIter] = force[qmmsg->coord[i].id];
02597             forceIter++;
02598         }
02599         
02600         delete qmmsg;
02601         
02602         if (pntChrgCoordMsgs != NULL) {
02603             DebugM(1,"Loading PntChrg forces.\n");
02604             
02605             for ( int i=0; i < pcmsg->numAtoms; ++i ) {
02606                 fmsg->force[forceIter] = force[pcmsg->coord[i].id];
02607                 forceIter++;
02608             }
02609             
02610             delete pcmsg;
02611         }
02612         
02613         DebugM(1,"A total of " << forceIter << " forces were loaded." << std::endl);
02614         
02615         for ( int i=0; i < subsArray.size(); ++i ) {
02616             fmsg->lssDat[i] = subsArray[i];
02617         }
02618         
02619         #ifdef DEBUG_QM
02620         if (subsArray.size() > 0)
02621             DebugM(3,"A total of " << subsArray.size() << " LSS data structures were loaded." << std::endl);
02622         #endif
02623         
02624         if ( ! j ) {
02625             fmsg->energy = totalEnergy;
02626             for ( int k=0; k<3; ++k )
02627                 for ( int l=0; l<3; ++l )
02628                     fmsg->virial[k][l] = totVirial[k][l];
02629         } else {
02630             fmsg->energy = 0;
02631             for ( int k=0; k<3; ++k )
02632                 for ( int l=0; l<3; ++l )
02633                     fmsg->virial[k][l] = 0;
02634         }
02635         
02636         DebugM(4,"Sending forces...\n");
02637         
02638         QMProxy[sourceNode].recvForce(fmsg);
02639         
02640     }
02641     
02642     DebugM(4,"All forces sent from node zero.\n");
02643 }

void ComputeQMMgr::recvForce ( QMForceMsg  ) 

Definition at line 2645 of file ComputeQM.C.

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

02645                                              {
02646     
02647     if (CkMyPe()) {
02648         for (int i=0; i<fmsg->numLssDat; i++) {
02649             subsArray.add(fmsg->lssDat[i]) ;
02650         }
02651     }
02652     
02653     QMCompute->saveResults(fmsg);
02654 }

void ComputeQMMgr::recvFullQM ( QMCoordMsg  ) 

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  ) 

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, 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, Node::Object(), PatchMap::Object(), open_dcd_write(), SimParameters::outputFilename, PCMODECUSTOMSEL, PCMODEUPDATEPOS, PCMODEUPDATESEL, 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  ) 

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, QMPntChrgMsg::coord, COULOMB, custom_ComputeQMAtom_Less(), SimParameters::cutoff, QMGrpCalcMsg::cutoff, QMGrpCalcMsg::data, StringList::data, DebugM, SimParameters::dielectric, QMAtomData::dist, QMAtomData::element, endi(), QMGrpCalcMsg::execPath, f, 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(), if(), 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(), 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  ) 

Definition at line 2422 of file ComputeQM.C.

References QMForce::charge, DebugM, endi(), QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, f, 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 Thu Jun 21 01:17:19 2018 for NAMD by  doxygen 1.4.7