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 400 of file ComputeQM.C.


Constructor & Destructor Documentation

ComputeQMMgr::ComputeQMMgr (  ) 

Definition at line 547 of file ComputeQM.C.

00547                            :
00548   QMProxy(thisgroup), QMCompute(0), numSources(0), numArrivedQMMsg(0), 
00549   numArrivedPntChrgMsg(0), numRecQMRes(0), qmCoordMsgs(0), pntChrgCoordMsgs(0),
00550   qmCoord(0), force(0), numAtoms(0), dcdOutFile(0), outputData(0), pcIDListSize(0) {
00551       
00552   CkpvAccess(BOCclass_group).computeQMMgr = thisgroup;
00553   
00554 }

ComputeQMMgr::~ComputeQMMgr (  ) 

Definition at line 556 of file ComputeQM.C.

References close_dcd_write().

00556                             {
00557     
00558     if (qmCoordMsgs != 0) {
00559         for ( int i=0; i<numSources; ++i ) {
00560             if (qmCoordMsgs[i] != 0)
00561                 delete qmCoordMsgs[i];
00562         }
00563         delete [] qmCoordMsgs;
00564     }
00565     
00566     if (pntChrgCoordMsgs != 0) {
00567         for ( int i=0; i<numSources; ++i ) {
00568             if (pntChrgCoordMsgs[i] != 0)
00569                 delete pntChrgCoordMsgs[i];
00570         }
00571         delete [] pntChrgCoordMsgs;
00572     }
00573     
00574     if (qmCoord != 0)
00575         delete [] qmCoord;
00576     
00577     if (force != 0)
00578         delete [] force;
00579     
00580     
00581     if (dcdOutFile != 0)
00582         close_dcd_write(dcdOutFile);
00583     
00584     if (dcdPosOutFile != 0)
00585         close_dcd_write(dcdPosOutFile);
00586     
00587     if (outputData != 0)
00588         delete [] outputData;
00589     
00590     if (lssPos != NULL)
00591         delete [] lssPos;
00592 }


Member Function Documentation

void ComputeQMMgr::calcMOPAC ( QMGrpCalcMsg  ) 

Definition at line 2704 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(), 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.

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

void ComputeQMMgr::calcORCA ( QMGrpCalcMsg  ) 

Definition at line 3531 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMForce::charge, QMAtomData::charge, QMGrpCalcMsg::charge, QMGrpCalcMsg::configLines, QMGrpCalcMsg::constants, QMGrpCalcMsg::data, DebugM, QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, QMGrpCalcMsg::execPath, QMGrpResMsg::force, QMGrpCalcMsg::grpIndx, QMGrpResMsg::grpIndx, QMAtomData::id, iERROR(), if(), iout, j, Vector::length(), QMGrpCalcMsg::multiplicity, NAMD_die(), 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, QMAtomData::type, QMGrpResMsg::virial, Vector::x, x, Vector::y, y, Vector::z, and z.

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

void ComputeQMMgr::calcUSR ( QMGrpCalcMsg  ) 

Definition at line 4362 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMForce::charge, QMAtomData::charge, QMGrpCalcMsg::constants, QMGrpCalcMsg::data, DebugM, QMGrpResMsg::energyCorr, QMGrpResMsg::energyOrig, QMGrpCalcMsg::execPath, QMGrpResMsg::force, QMGrpCalcMsg::grpIndx, QMGrpResMsg::grpIndx, QMAtomData::id, iERROR(), if(), iout, j, Vector::length(), NAMD_die(), 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, QMAtomData::type, QMGrpResMsg::virial, Vector::x, x, Vector::y, y, Vector::z, and z.

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

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

Definition at line 427 of file ComputeQM.C.

Referenced by lssSubs().

00427 { return subsArray; } ;

void ComputeQMMgr::procQMRes (  ) 

Definition at line 2406 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.

02406                              {
02407     
02408     // Writes a DCD file with the charges of all QM atoms at a frequency 
02409     // defined by the user in qmOutFreq.
02410     if ( simParams->qmOutFreq > 0 && 
02411          timeStep % simParams->qmOutFreq == 0 ) {
02412         
02413         iout << iINFO << "Writing QM charge output at step " 
02414         << timeStep <<  "\n" << endi;
02415     
02416         Real *x = outputData, 
02417              *y = outputData + molPtr->get_numQMAtoms(), 
02418              *z = outputData + 2*molPtr->get_numQMAtoms();
02419         
02420         for (int qmIt=0; qmIt<numQMAtms; qmIt++){
02421             x[qmIt] = qmCoord[qmIt].id;
02422             y[qmIt] = force[qmCoord[qmIt].id].charge ;
02423             z[qmIt] = 0;
02424         }
02425         
02426         write_dcdstep(dcdOutFile, numQMAtms, x, y, z, 0) ;
02427     }
02428     
02429     // Writes a DCD file with the positions of all QM atoms at a frequency 
02430     // defined by the user in qmPosOutFreq.
02431     if ( simParams->qmPosOutFreq > 0 && 
02432          timeStep % simParams->qmPosOutFreq == 0 ) {
02433         
02434         iout << iINFO << "Writing QM position output at step " 
02435         << timeStep <<  "\n" << endi;
02436         
02437         SortedArray<idIndxStr> idIndx;
02438         
02439         for(int i=0; i<numQMAtms;i++) {
02440             idIndx.insert( idIndxStr(qmCoord[i].id, i) );
02441         }
02442         idIndx.sort();
02443         
02444         Real *x = outputData, 
02445              *y = outputData + molPtr->get_numQMAtoms(), 
02446              *z = outputData + 2*molPtr->get_numQMAtoms();
02447         
02448         for (int qmIt=0; qmIt<numQMAtms; qmIt++){
02449             x[qmIt] = qmCoord[idIndx[qmIt].indx].position.x;
02450             y[qmIt] = qmCoord[idIndx[qmIt].indx].position.y;
02451             z[qmIt] = qmCoord[idIndx[qmIt].indx].position.z;
02452         }
02453         
02454         write_dcdstep(dcdPosOutFile, numQMAtms, x, y, z, 0) ;
02455     }
02456     
02457     // distribute forces
02458     DebugM(4,"Distributing QM forces for all ranks.\n");
02459     for ( int j=0; j < numSources; ++j ) {
02460         
02461         DebugM(1,"Sending forces and charges to source " << j << std::endl);
02462         
02463         QMCoordMsg *qmmsg = 0;
02464         QMPntChrgMsg *pcmsg = 0 ;
02465         
02466         int totForces = 0;
02467         int sourceNode = -1;
02468         
02469         if (pntChrgCoordMsgs == NULL) {
02470             
02471             qmmsg = qmCoordMsgs[j];
02472             qmCoordMsgs[j] = 0;
02473             
02474             totForces = qmmsg->numAtoms ;
02475             
02476             sourceNode = qmmsg->sourceNode;
02477         }
02478         else {
02479             pcmsg = pntChrgCoordMsgs[j];
02480             pntChrgCoordMsgs[j] = 0;
02481             
02482             sourceNode = pcmsg->sourceNode;
02483             
02484             // Since we receive two different messages from nodes, there is no 
02485             // guarantee the two sets of messages will come in the same order.
02486             // Therefore, we match the messages by comaring their sourceNodes.
02487             for (int aux=0; aux<numSources; aux++) {
02488                 
02489                 if (qmCoordMsgs[aux] == 0)
02490                     continue;
02491                 
02492                 qmmsg = qmCoordMsgs[aux];
02493                 
02494                 if (qmmsg->sourceNode == sourceNode) {
02495                     qmCoordMsgs[aux] = 0;
02496                     break;
02497                 }
02498             }
02499             
02500             DebugM(1,"Building force mesage for rank " 
02501             << pcmsg->sourceNode << std::endl);
02502             
02503             totForces = qmmsg->numAtoms + pcmsg->numAtoms;
02504         }
02505         
02506         QMForceMsg *fmsg = new (totForces, subsArray.size(), 0) QMForceMsg;
02507         fmsg->PMEOn = simParams->PMEOn;
02508         fmsg->numForces = totForces;
02509         fmsg->numLssDat = subsArray.size();
02510         
02511         DebugM(1,"Loading QM forces.\n");
02512         
02513         // This iterator is used in BOTH loops!
02514         int forceIter = 0;
02515         
02516         for ( int i=0; i < qmmsg->numAtoms; ++i ) {
02517             fmsg->force[forceIter] = force[qmmsg->coord[i].id];
02518             forceIter++;
02519         }
02520         
02521         delete qmmsg;
02522         
02523         if (pntChrgCoordMsgs != NULL) {
02524             DebugM(1,"Loading PntChrg forces.\n");
02525             
02526             for ( int i=0; i < pcmsg->numAtoms; ++i ) {
02527                 fmsg->force[forceIter] = force[pcmsg->coord[i].id];
02528                 forceIter++;
02529             }
02530             
02531             delete pcmsg;
02532         }
02533         
02534         DebugM(1,"A total of " << forceIter << " forces were loaded." << std::endl);
02535         
02536         for ( int i=0; i < subsArray.size(); ++i ) {
02537             fmsg->lssDat[i] = subsArray[i];
02538         }
02539         
02540         #ifdef DEBUG_QM
02541         if (subsArray.size() > 0)
02542             DebugM(3,"A total of " << subsArray.size() << " LSS data structures were loaded." << std::endl);
02543         #endif
02544         
02545         if ( ! j ) {
02546             fmsg->energy = totalEnergy;
02547             for ( int k=0; k<3; ++k )
02548                 for ( int l=0; l<3; ++l )
02549                     fmsg->virial[k][l] = totVirial[k][l];
02550         } else {
02551             fmsg->energy = 0;
02552             for ( int k=0; k<3; ++k )
02553                 for ( int l=0; l<3; ++l )
02554                     fmsg->virial[k][l] = 0;
02555         }
02556         
02557         DebugM(4,"Sending forces...\n");
02558         
02559         QMProxy[sourceNode].recvForce(fmsg);
02560         
02561     }
02562     
02563     DebugM(4,"All forces sent from node zero.\n");
02564 }

void ComputeQMMgr::recvForce ( QMForceMsg  ) 

Definition at line 2566 of file ComputeQM.C.

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

02566                                              {
02567     
02568     if (CkMyPe()) {
02569         for (int i=0; i<fmsg->numLssDat; i++) {
02570             subsArray.add(fmsg->lssDat[i]) ;
02571         }
02572     }
02573     
02574     QMCompute->saveResults(fmsg);
02575 }

void ComputeQMMgr::recvFullQM ( QMCoordMsg  ) 

Definition at line 1280 of file ComputeQM.C.

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

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

void ComputeQMMgr::recvPartQM ( QMCoordMsg  ) 

Definition at line 817 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().

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

void ComputeQMMgr::recvPntChrg ( QMPntChrgMsg  ) 

Definition at line 1832 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, 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, 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, QMAtomData::type, and Vector::unit().

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

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

Definition at line 405 of file ComputeQM.C.

00405 { QMCompute = c; }

void ComputeQMMgr::storeQMRes ( QMGrpResMsg  ) 

Definition at line 2343 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.

02343                                                  {
02344     
02345 //     iout << iINFO << "Storing QM results for region " << resMsg->grpIndx  
02346 //     << " (ID: "  << grpID[resMsg->grpIndx] 
02347 //     << ") with original energy: " << endi;
02348 //     std::cout << std::fixed << std::setprecision(6) << resMsg->energyOrig << endi;
02349 //     iout << " Kcal/mol\n" << endi;
02350     
02351 //     if (resMsg->energyCorr != resMsg->energyOrig) {
02352 //         iout << iINFO << "PME corrected energy: " << endi;
02353 //         std::cout << std::fixed << std::setprecision(6) << resMsg->energyCorr << endi;
02354 //         iout << " Kcal/mol\n" << endi;
02355 //     }
02356     
02357     if ( (timeStep % simParams->qmEnergyOutFreq ) == 0 ) {
02358         iout << QMETITLE(timeStep);
02359         iout << FORMAT(grpID[resMsg->grpIndx]);
02360         iout << FORMAT(resMsg->energyOrig);
02361         if (resMsg->energyCorr != resMsg->energyOrig) iout << FORMAT(resMsg->energyCorr);
02362         iout << "\n\n" << endi;
02363     }
02364     
02365     totalEnergy += resMsg->energyCorr ;
02366     
02367     for ( int k=0; k<3; ++k )
02368         for ( int l=0; l<3; ++l )
02369             totVirial[k][l] += resMsg->virial[k][l];
02370     
02371     QMForce *fres = resMsg->force ;
02372     Real qmTotalCharge = 0;
02373     
02374     for (int i=0; i<resMsg->numForces; i++) {
02375         
02376         force[ fres[i].id ].force += fres[i].force;
02377         
02378         // Indicates the result is a QM atom, and we should get its updated charge.
02379         if (fres[i].replace == 1) {
02380             force[ fres[i].id ].charge =  fres[i].charge;
02381             qmTotalCharge += fres[i].charge;
02382         }
02383     }
02384     
02385     if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
02386         qmTotalCharge = roundf(qmTotalCharge) ;
02387     }
02388     
02389     // In case we are calculating cSMD forces, apply them now.
02390     if (cSMDon) {
02391         for ( int i=0; i < cSMDindxLen[resMsg->grpIndx]; i++ ) {
02392             int cSMDit = cSMDindex[resMsg->grpIndx][i];
02393             force[ cSMDpairs[cSMDit][0] ].force += cSMDForces[cSMDit];
02394         }
02395     }
02396 
02397     DebugM(4,"QM total charge received is " << qmTotalCharge << std::endl);
02398     
02399     DebugM(4,"Current accumulated energy is " << totalEnergy << std::endl);
02400     
02401     numRecQMRes++;
02402     
02403     delete resMsg;
02404 }


The documentation for this class was generated from the following file:
Generated on Tue Nov 21 01:17:18 2017 for NAMD by  doxygen 1.4.7