version 1.3 | version 1.4 |
---|
| |
class QMGrpResMsg : public CMessage_QMGrpResMsg { | class QMGrpResMsg : public CMessage_QMGrpResMsg { |
public: | public: |
int grpIndx; | int grpIndx; |
BigReal energy; | BigReal energyOrig; // Original QM Energy |
| BigReal energyCorr; // Corrected energy due to PME |
BigReal virial[3][3]; | BigReal virial[3][3]; |
int numForces; | int numForces; |
QMForce *force; | QMForce *force; |
| |
| |
iout << iINFO << "Storing QM results for region " << resMsg->grpIndx | iout << iINFO << "Storing QM results for region " << resMsg->grpIndx |
<< " (ID: " << grpID[resMsg->grpIndx] | << " (ID: " << grpID[resMsg->grpIndx] |
<< ") with energy: " << endi; | << ") with original energy: " << endi; |
std::cout << std::fixed << std::setprecision(6) << resMsg->energy << endi; | std::cout << std::fixed << std::setprecision(6) << resMsg->energyOrig << endi; |
iout << " Kcal/mol\n" << endi; | iout << " Kcal/mol\n" << endi; |
| |
totalEnergy += resMsg->energy ; | if (resMsg->energyCorr != resMsg->energyOrig) { |
| iout << iINFO << "PME corrected energy: " << endi; |
| std::cout << std::fixed << std::setprecision(6) << resMsg->energyCorr << endi; |
| iout << " Kcal/mol\n" << endi; |
| } |
| |
| totalEnergy += resMsg->energyCorr ; |
| |
for ( int k=0; k<3; ++k ) | for ( int k=0; k<3; ++k ) |
for ( int l=0; l<3; ++l ) | for ( int l=0; l<3; ++l ) |
| |
| |
} | } |
| |
DebugM(1,"Placing forces on force boxes.\n"); | DebugM(1,"Placing forces on force boxes in rank " |
| << CkMyPe() << std::endl); |
| |
// Places the received forces on the force array for each patch. | // Places the received forces on the force array for each patch. |
int homeIndxIter = 0; | int homeIndxIter = 0; |
| |
| |
} | } |
| |
| DebugM(1,"Forces placed on force boxes in rank " |
| << CkMyPe() << std::endl); |
| |
if (fmsg->PMEOn) { | if (fmsg->PMEOn) { |
| |
| DebugM(1,"PME ON! Accessing PMEmgr in rank " << CkMyPe() << std::endl); |
| |
ComputePmeMgr *mgrP = CProxy_ComputePmeMgr::ckLocalBranch( | ComputePmeMgr *mgrP = CProxy_ComputePmeMgr::ckLocalBranch( |
CkpvAccess(BOCclass_group).computePmeMgr) ; | CkpvAccess(BOCclass_group).computePmeMgr) ; |
| |
| |
QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg; | QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg; |
resMsg->grpIndx = msg->grpIndx; | resMsg->grpIndx = msg->grpIndx; |
resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs; | resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs; |
resMsg->energy = 0; | resMsg->energyOrig = 0; |
| resMsg->energyCorr = 0; |
for ( int k=0; k<3; ++k ) | for ( int k=0; k<3; ++k ) |
for ( int l=0; l<3; ++l ) | for ( int l=0; l<3; ++l ) |
resMsg->virial[k][l] = 0; | resMsg->virial[k][l] = 0; |
| |
| |
// We have to convert the notation from FORTRAN double precision to | // We have to convert the notation from FORTRAN double precision to |
// the natural, normal, reasonable and not terribly out dated format. | // the natural, normal, reasonable and not terribly out dated format. |
resMsg->energy = strtod(strEnergy, &endPtr); | resMsg->energyOrig = strtod(strEnergy, &endPtr); |
if ( *endPtr == 'D' ) { | if ( *endPtr == 'D' ) { |
*endPtr = 'e'; | *endPtr = 'e'; |
resMsg->energy = strtod(strEnergy, &endPtr); | resMsg->energyOrig = strtod(strEnergy, &endPtr); |
} | } |
| |
// In MOPAC, the total energy is given in EV, so we convert to Kcal/mol | // In MOPAC, the total energy is given in EV, so we convert to Kcal/mol |
resMsg->energy *= 23.060 ; | resMsg->energyOrig *= 23.060 ; |
| |
// DebugM(4,"Reading QM energy from file: " << resMsg->energy << "\n"); | // DebugM(4,"Reading QM energy from file: " << resMsg->energyOrig << "\n"); |
| |
| resMsg->energyCorr = resMsg->energyOrig; |
| |
continue; | continue; |
} | } |
| |
resForce[j].force -= -1*fixForce ; | resForce[j].force -= -1*fixForce ; |
| |
// The energy is *subtracted* from the total energy calculated here. | // The energy is *subtracted* from the total energy calculated here. |
resMsg->energy -= energy; | resMsg->energyCorr -= energy; |
} | } |
} | } |
| |
// DebugM(4,"Corrected energy QM-QM interactions: " << resMsg->energy << "\n"); | // DebugM(4,"Corrected energy QM-QM interactions: " << resMsg->energyCorr << "\n"); |
| |
pcIndx = msg->numQMAtoms; | pcIndx = msg->numQMAtoms; |
// We only loop over point charges from real atoms, ignoring the ones | // We only loop over point charges from real atoms, ignoring the ones |
| |
fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ; | fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ; |
| |
// The energy is *subtracted* from the total energy calculated here. | // The energy is *subtracted* from the total energy calculated here. |
resMsg->energy -= energy; | resMsg->energyCorr -= energy; |
} | } |
| |
// The force is *subtracted* from the total force acting on | // The force is *subtracted* from the total force acting on |
| |
QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg; | QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg; |
resMsg->grpIndx = msg->grpIndx; | resMsg->grpIndx = msg->grpIndx; |
resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs; | resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs; |
resMsg->energy = 0; | resMsg->energyOrig = 0; |
| resMsg->energyCorr = 0; |
for ( int k=0; k<3; ++k ) | for ( int k=0; k<3; ++k ) |
for ( int l=0; l<3; ++l ) | for ( int l=0; l<3; ++l ) |
resMsg->virial[k][l] = 0; | resMsg->virial[k][l] = 0; |
| |
fgets(line, lineLen, outputFile); | fgets(line, lineLen, outputFile); |
} | } |
| |
iret = fscanf(outputFile,"%lf\n", &resMsg->energy); | iret = fscanf(outputFile,"%lf\n", &resMsg->energyOrig); |
if ( iret != 1 ) { | if ( iret != 1 ) { |
NAMD_die("Error reading QM forces file."); | NAMD_die("Error reading QM forces file."); |
} | } |
// iout << "Energy in step (Hartree): " << resMsg->energy << "\n" << endi; | // iout << "Energy in step (Hartree): " << resMsg->energyOrig << "\n" << endi; |
// All energies are given in Eh (Hartree) | // All energies are given in Eh (Hartree) |
// NAMD needs energies in kcal/mol | // NAMD needs energies in kcal/mol |
// The conversion factor is 627.509469 | // The conversion factor is 627.509469 |
resMsg->energy *= 627.509469; | resMsg->energyOrig *= 627.509469; |
// iout << "Energy in step (Kcal/mol): " << resMsg->energy << "\n" << endi; | // iout << "Energy in step (Kcal/mol): " << resMsg->energyOrig << "\n" << endi; |
| |
| resMsg->energyCorr = resMsg->energyOrig; |
| |
// skip lines before gradient | // skip lines before gradient |
for (int i = 0; i < 3; i++) { | for (int i = 0; i < 3; i++) { |
| |
resForce[j].force -= -1*fixForce; | resForce[j].force -= -1*fixForce; |
| |
// The energy is *subtracted* from the total energy calculated here. | // The energy is *subtracted* from the total energy calculated here. |
resMsg->energy -= energy; | resMsg->energyCorr -= energy; |
} | } |
} | } |
| |
| |
fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ; | fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ; |
| |
// The energy is *subtracted* from the total energy calculated here. | // The energy is *subtracted* from the total energy calculated here. |
resMsg->energy -= energy; | resMsg->energyCorr -= energy; |
| |
} | } |
| |
| |
QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg; | QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg; |
resMsg->grpIndx = msg->grpIndx; | resMsg->grpIndx = msg->grpIndx; |
resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs; | resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs; |
resMsg->energy = 0; | resMsg->energyOrig = 0; |
| resMsg->energyCorr = 0; |
for ( int k=0; k<3; ++k ) | for ( int k=0; k<3; ++k ) |
for ( int l=0; l<3; ++l ) | for ( int l=0; l<3; ++l ) |
resMsg->virial[k][l] = 0; | resMsg->virial[k][l] = 0; |
| |
// Reads the data form the output file created by the QM software. | // Reads the data form the output file created by the QM software. |
// Gradients over the QM atoms, and Charges for QM atoms will be read. | // Gradients over the QM atoms, and Charges for QM atoms will be read. |
| |
iret = fscanf(outputFile,"%lf\n", &resMsg->energy); | iret = fscanf(outputFile,"%lf\n", &resMsg->energyOrig); |
if ( iret != 1 ) { | if ( iret != 1 ) { |
NAMD_die("Error reading QM forces file."); | NAMD_die("Error reading energy from QM results file."); |
} | } |
| |
| resMsg->energyCorr = resMsg->energyOrig; |
| |
size_t atmIndx; | size_t atmIndx; |
double localForce[3]; | double localForce[3]; |
double localCharge; | double localCharge; |
for (atmIndx = 0; atmIndx < msg->numAllAtoms; atmIndx++) { | for (atmIndx = 0; atmIndx < msg->numAllAtoms; atmIndx++) { |
| |
| |
| |
iret = fscanf(outputFile,"%lf %lf %lf %lf\n", | iret = fscanf(outputFile,"%lf %lf %lf %lf\n", |
localForce+0, | localForce+0, |
localForce+1, | localForce+1, |
localForce+2, | localForce+2, |
&localCharge); | &localCharge); |
if ( iret != 1 ) { | if ( iret != 4 ) { |
NAMD_die("Error reading QM forces file."); | NAMD_die("Error reading gradient and charge from QM results file."); |
} | } |
| |
// If we are reading charges and forces on QM atoms, store | // If we are reading charges and forces on QM atoms, store |
| |
resForce[j].force -= -1*fixForce; | resForce[j].force -= -1*fixForce; |
| |
// The energy is *subtracted* from the total energy calculated here. | // The energy is *subtracted* from the total energy calculated here. |
resMsg->energy -= energy; | resMsg->energyCorr -= energy; |
} | } |
} | } |
| |
| |
fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ; | fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ; |
| |
// The energy is *subtracted* from the total energy calculated here. | // The energy is *subtracted* from the total energy calculated here. |
resMsg->energy -= energy; | resMsg->energyCorr -= energy; |
| |
} | } |
| |