Difference for src/ComputeQM.C from version 1.3 to 1.4

version 1.3version 1.4
Line 176
Line 176
 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;
Line 2234
Line 2235
          
     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 )
Line 2508
Line 2515
                  
     }     }
          
     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;
Line 2532
Line 2540
                  
     }     }
          
      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) ;
                  
Line 2771
Line 2785
     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;
Line 2827
Line 2842
                          
             // 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;
         }         }
Line 3157
Line 3174
                 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 
Line 3204
Line 3221
                 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
Line 3476
Line 3493
     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;
Line 3535
Line 3553
         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++) {
Line 3936
Line 3956
                 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;
             }             }
         }         }
                  
Line 3981
Line 4001
                 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;
                                  
             }             }
                          
Line 4191
Line 4211
     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;
Line 4225
Line 4246
     // 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
Line 4454
Line 4475
                 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;
             }             }
         }         }
                  
Line 4499
Line 4520
                 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;
                                  
             }             }
                          

Legend:
Removed in v.1.3 
changed lines
 Added in v.1.4


Made by using version 1.53 of cvs2html