NAMD
Public Member Functions | List of all members
ComputeQMMgr Class Reference
Inheritance diagram for ComputeQMMgr:

Public Member Functions

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

Detailed Description

Definition at line 401 of file ComputeQM.C.

Constructor & Destructor Documentation

ComputeQMMgr::ComputeQMMgr ( )

Definition at line 548 of file ComputeQM.C.

548  :
549  QMProxy(thisgroup), QMCompute(0), numSources(0), numArrivedQMMsg(0),
550  numArrivedPntChrgMsg(0), numRecQMRes(0), qmCoordMsgs(0), pntChrgCoordMsgs(0),
551  qmCoord(0), force(0), numAtoms(0), dcdOutFile(0), outputData(0), pcIDListSize(0) {
552 
553  CkpvAccess(BOCclass_group).computeQMMgr = thisgroup;
554 
555 }
ComputeQMMgr::~ComputeQMMgr ( )

Definition at line 557 of file ComputeQM.C.

References close_dcd_write().

557  {
558 
559  if (qmCoordMsgs != 0) {
560  for ( int i=0; i<numSources; ++i ) {
561  if (qmCoordMsgs[i] != 0)
562  delete qmCoordMsgs[i];
563  }
564  delete [] qmCoordMsgs;
565  }
566 
567  if (pntChrgCoordMsgs != 0) {
568  for ( int i=0; i<numSources; ++i ) {
569  if (pntChrgCoordMsgs[i] != 0)
570  delete pntChrgCoordMsgs[i];
571  }
572  delete [] pntChrgCoordMsgs;
573  }
574 
575  if (qmCoord != 0)
576  delete [] qmCoord;
577 
578  if (force != 0)
579  delete [] force;
580 
581 
582  if (dcdOutFile != 0)
583  close_dcd_write(dcdOutFile);
584 
585  if (dcdPosOutFile != 0)
586  close_dcd_write(dcdPosOutFile);
587 
588  if (outputData != 0)
589  delete [] outputData;
590 
591  if (lssPos != NULL)
592  delete [] lssPos;
593 }
close_dcd_write
void close_dcd_write(int fd)
Definition: dcdlib.C:1063

Member Function Documentation

void ComputeQMMgr::calcMOPAC ( QMGrpCalcMsg *  msg)

Definition at line 2801 of file ComputeQM.C.

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

2802 {
2803  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
2804  FILE *inputFile,*outputFile,*chrgFile;
2805  int iret;
2806 
2807  const size_t lineLen = 256;
2808  char *line = new char[lineLen];
2809 
2810  std::string qmCommand, inputFileName, outputFileName, pntChrgFileName;
2811 
2812  // For coulomb calculation
2813  BigReal constants = msg->constants;
2814 
2815  double gradient[3];
2816 
2817  DebugM(4,"Running MOPAC on PE " << CkMyPe() << std::endl);
2818 
2819  QMAtomData *atmP = msg->data ;
2820  QMAtomData *pcP = msg->data + msg->numAllAtoms ;
2821 
2822  // We create a pointer for PME correction, which may point to
2823  // a copy of the original point charge array, with unchanged charges,
2824  // or points to the original array in case no switching or charge
2825  // scheme is being used.
2826  QMAtomData *pcPpme = NULL;
2827  if (msg->switching) {
2828 
2829  if (msg->PMEOn)
2830  pcPpme = new QMAtomData[msg->numRealPntChrgs];
2831 
2832  pntChrgSwitching(msg, pcPpme) ;
2833  } else {
2834  pcPpme = pcP;
2835  }
2836 
2837  int retVal = 0;
2838  struct stat info;
2839 
2840  // For each QM group, create a subdirectory where files will be palced.
2841  std::string baseDir(msg->baseDir);
2842  std::ostringstream itosConv ;
2843  if ( CmiNumPartitions() > 1 ) {
2844  baseDir.append("/") ;
2845  itosConv << CmiMyPartition() ;
2846  baseDir += itosConv.str() ;
2847  itosConv.str("");
2848  itosConv.clear() ;
2849 
2850  if (stat(msg->baseDir, &info) != 0 ) {
2851  CkPrintf( "Node %d cannot access directory %s\n",
2852  CkMyPe(), baseDir.c_str() );
2853  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
2854  }
2855  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
2856  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
2857  retVal = mkdir(baseDir.c_str(), S_IRWXU);
2858  }
2859 
2860  }
2861  baseDir.append("/") ;
2862  itosConv << msg->grpIndx ;
2863  baseDir += itosConv.str() ;
2864 
2865  if (stat(msg->baseDir, &info) != 0 ) {
2866  CkPrintf( "Node %d cannot access directory %s\n",
2867  CkMyPe(), baseDir.c_str() );
2868  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
2869  }
2870  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
2871  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
2872  retVal = mkdir(baseDir.c_str(), S_IRWXU);
2873  }
2874 
2875  #ifdef DEBUG_QM
2876  Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
2877  #endif
2878 
2879  inputFileName.clear();
2880  inputFileName.append(baseDir.c_str()) ;
2881  inputFileName.append("/qmmm_") ;
2882  inputFileName += itosConv.str() ;
2883  inputFileName.append(".input") ;
2884 
2885  // Opens file for coordinate and parameter input
2886  inputFile = fopen(inputFileName.c_str(),"w");
2887  if ( ! inputFile ) {
2888  iout << iERROR << "Could not open input file for writing: "
2889  << inputFileName << "\n" << endi ;
2890  NAMD_err("Error writing QM input file.");
2891  }
2892 
2893  // Builds the command that will be executed
2894  qmCommand.clear();
2895  qmCommand.append("cd ");
2896  qmCommand.append(baseDir);
2897  qmCommand.append(" ; ");
2898  qmCommand.append(msg->execPath) ;
2899  qmCommand.append(" ") ;
2900  qmCommand.append(inputFileName) ;
2901  qmCommand.append(" > /dev/null 2>&1") ;
2902 
2903  // Builds the file name where MOPAC will place the gradient
2904  // This is also relative to the input file name.
2905  outputFileName = inputFileName ;
2906  outputFileName.append(".aux") ;
2907 
2908  if (msg->numAllPntChrgs) {
2909  // Builds the file name where we will write the point charges.
2910  pntChrgFileName.clear();
2911  pntChrgFileName.append(baseDir.c_str()) ;
2912  pntChrgFileName.append("/mol.in") ;
2913 
2914  chrgFile = fopen(pntChrgFileName.c_str(),"w");
2915  if ( ! chrgFile ) {
2916  iout << iERROR << "Could not open charge file for writing: "
2917  << pntChrgFileName << "\n" << endi ;
2918  NAMD_err("Error writing point charge file.");
2919  }
2920 
2921  iret = fprintf(chrgFile,"\n%d %d\n",
2922  msg->numQMAtoms, msg->numAllAtoms - msg->numQMAtoms);
2923  if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
2924  }
2925 
2926  // writes configuration lines to the input file.
2927  std::stringstream ss(msg->configLines) ;
2928  std::string confLineStr;
2929  while (std::getline(ss, confLineStr) ) {
2930  confLineStr.append("\n");
2931  iret = fprintf(inputFile,confLineStr.c_str());
2932  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
2933  }
2934 
2935 
2936  iret = fprintf(inputFile,"\n");
2937  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
2938 
2939  DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file "
2940  << inputFileName.c_str() << " and " << msg->numAllPntChrgs
2941  << " point charges in file " << pntChrgFileName.c_str() << "\n");
2942 
2943  // write QM and dummy atom coordinates to input file and
2944  // MM electric field from MM point charges.
2945  for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
2946 
2947  double x = atmP->position.x;
2948  double y = atmP->position.y;
2949  double z = atmP->position.z;
2950 
2951  iret = fprintf(inputFile,"%s %f 1 %f 1 %f 1\n",
2952  atmP->element,x,y,z);
2953  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
2954 
2955  if (msg->numAllPntChrgs) {
2956  BigReal phi = 0;
2957 
2958  // The Electrostatic Potential is calculated according to
2959  // the QMMM Keyword documentation for MOPAC
2960  // http://openmopac.net/manual/QMMM.html
2961  pcP = msg->data + msg->numAllAtoms ;
2962  for ( size_t j=0; j < msg->numAllPntChrgs; ++j, ++pcP ) {
2963 
2964  // In case of QM-MM bonds, the charge of the MM1 atom is ignored
2965  if (pcP->type == QMPCTYPE_IGNORE)
2966  continue;
2967 
2968  double charge = pcP->charge;
2969 
2970  double xMM = pcP->position.x;
2971  double yMM = pcP->position.y;
2972  double zMM = pcP->position.z;
2973 
2974  BigReal rij = sqrt((x-xMM)*(x-xMM) +
2975  (y-yMM)*(y-yMM) +
2976  (z-zMM)*(z-zMM) ) ;
2977 
2978  phi += charge/rij ;
2979  }
2980 
2981  // We use the same Coulomb constant used in the rest of NAMD
2982  // instead of the suggested rounded 332 suggested by MOPAC.
2983  phi = phi*constants ;
2984 
2985  iret = fprintf(chrgFile,"%s %f %f %f %f\n",
2986  atmP->element,x,y,z, phi);
2987  if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
2988  }
2989  }
2990 
2991  DebugM(4,"Closing input and charge file\n");
2992 
2993  if (msg->numAllPntChrgs)
2994  fclose(chrgFile);
2995 
2996  fclose(inputFile);
2997 
2998  if (msg->prepProcOn) {
2999 
3000  std::string prepProc(msg->prepProc) ;
3001  prepProc.append(" ") ;
3002  prepProc.append(inputFileName) ;
3003  iret = system(prepProc.c_str());
3004  if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
3005  if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
3006  }
3007 
3008  // runs QM command
3009  DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
3010  iret = system(qmCommand.c_str());
3011 
3012  if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
3013  if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
3014 
3015  if (msg->secProcOn) {
3016 
3017  std::string secProc(msg->secProc) ;
3018  secProc.append(" ") ;
3019  secProc.append(inputFileName) ;
3020  itosConv.str("");
3021  itosConv.clear() ;
3022  itosConv << msg->timestep ;
3023  secProc.append(" ") ;
3024  secProc += itosConv.str() ;
3025 
3026  iret = system(secProc.c_str());
3027  if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
3028  if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
3029  }
3030 
3031  // remove coordinate file
3032 // iret = remove(inputFileName);
3033 // if ( iret ) { NAMD_err(0); }
3034 
3035  // remove coordinate file
3036 // iret = remove(pntChrgFileName);
3037 // if ( iret ) { NAMD_err(0); }
3038 
3039  // opens output file
3040  DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
3041  outputFile = fopen(outputFileName.c_str(),"r");
3042  if ( ! outputFile ) {
3043  iout << iERROR << "Could not find QM output file!\n" << endi;
3044  NAMD_err(0);
3045  }
3046 
3047  // Resets the pointers.
3048  atmP = msg->data ;
3049  pcP = msg->data + msg->numAllAtoms ;
3050 
3051  // Allocates the return message.
3052  QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
3053  resMsg->grpIndx = msg->grpIndx;
3054  resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
3055  resMsg->energyOrig = 0;
3056  resMsg->energyCorr = 0;
3057  for ( int k=0; k<3; ++k )
3058  for ( int l=0; l<3; ++l )
3059  resMsg->virial[k][l] = 0;
3060  QMForce *resForce = resMsg->force ;
3061 
3062  // We split the data into two pointers so we can "skip" the dummy atoms
3063  // which have no representation in NAMD.
3064  for (int i=0; i<resMsg->numForces; i++) {
3065  resForce[i].force = 0;
3066  resForce[i].charge = 0 ;
3067  if (i < msg->numQMAtoms) {
3068  // We use the replace field to indicate QM atoms,
3069  // which will have charge information.
3070  resForce[i].replace = 1 ;
3071  resForce[i].id = atmP->id;
3072  atmP++;
3073  }
3074  else {
3075  // We use the replace field to indicate QM atoms,
3076  // which will have charge information.
3077  resForce[i].replace = 0 ;
3078  resForce[i].id = pcP->id;
3079  pcP++;
3080  }
3081  }
3082 
3083  // Resets the pointers.
3084  atmP = msg->data ;
3085  pcP = msg->data + msg->numAllAtoms ;
3086 
3087  // Reads the data form the output file created by the QM software.
3088  // Gradients over the QM atoms, and Charges for QM atoms will be read.
3089 
3090  size_t atmIndx = 0, gradCount = 0;
3091  Bool gradFields = false, chargeFields = false;
3092  Bool chargesRead = false, gradsRead = false;
3093  while ( fgets(line, lineLen, outputFile) != NULL){
3094  // We loop over the lines of the output file untill we find
3095  // the line that initiates the "atom charges" lines. Then
3096  // we read all lines and expect to get one or more charges
3097  // per line, separated by space(s), untill we find a line that
3098  // initiates the "gradients", then once more, we expect several
3099  // numbers separated by space(s). When the "overlap matrix"
3100  // string is found, we break the loop and stop reading the file.
3101 
3102  // if ( strstr(line,"TOTAL_ENERGY") != NULL ) {
3103  if ( strstr(line,"HEAT_OF_FORMATION") != NULL ) {
3104 
3105  char strEnergy[14], *endPtr ;
3106 
3107  //strncpy(strEnergy, line + 17, 13) ;
3108  strncpy(strEnergy, line + 28, 13) ;
3109  strEnergy[13] = '\0';
3110 
3111  // We have to convert the notation from FORTRAN double precision to
3112  // the natural, normal, reasonable and not terribly out dated format.
3113  resMsg->energyOrig = strtod(strEnergy, &endPtr);
3114  if ( *endPtr == 'D' ) {
3115  *endPtr = 'e';
3116  resMsg->energyOrig = strtod(strEnergy, &endPtr);
3117  }
3118 
3119  // In MOPAC, the total energy is given in EV, so we convert to Kcal/mol
3120 // resMsg->energyOrig *= 23.060 ; // We now read Heat of Formation, which is given in Kcal/Mol
3121 
3122 // DebugM(4,"Reading QM energy from file: " << resMsg->energyOrig << "\n");
3123 
3124  resMsg->energyCorr = resMsg->energyOrig;
3125 
3126  continue;
3127  }
3128 
3129  if ( strstr(line,"ATOM_CHARGES") != NULL ) {
3130  gradFields = false;
3131  chargeFields = true;
3132  atmIndx = 0;
3133 
3134  // If the user asked for charges NOT to be read from MOPAC,
3135  // we skip the charge reading step.
3136  if (msg->qmAtmChrgMode == QMCHRGNONE) {
3137  chargeFields = false;
3138  atmIndx = msg->numAllAtoms;
3139  }
3140  else {
3141  chargeFields = true;
3142  atmIndx = 0;
3143  }
3144 
3145  // Now we expect the following line(s) to have atom charges
3146  continue;
3147  }
3148 
3149  if ( strstr(line,"GRADIENTS") != NULL ) {
3150 
3151  // Now that all charges have been read, checks if the
3152  // numbers match
3153  if (atmIndx != msg->numAllAtoms) {
3154  NAMD_die("Error reading QM forces file. Wrong number of atom charges");
3155  }
3156 
3157  chargesRead = true;
3158 
3159  // Now we expect the following line(s) to have gradients
3160  chargeFields = false ;
3161  gradFields = true;
3162  gradCount = 0;
3163  atmIndx = 0;
3164 
3165  continue;
3166  }
3167 
3168  if ( strstr(line,"OVERLAP_MATRIX") != NULL ) {
3169 
3170  // Now that all gradients have been read, checks if the
3171  // numbers match
3172  if (atmIndx != msg->numAllAtoms) {
3173  NAMD_die("Error reading QM forces file. Wrong number of gradients");
3174  }
3175 
3176  gradsRead = true;
3177 
3178  // Nothing more to read from the ".aux" file
3179  break;
3180  }
3181 
3182  char result[20] ;
3183  size_t strIndx = 0;
3184 
3185  if (chargeFields) {
3186  while ((strIndx < (strlen(line)-9)) && (strlen(line)-1 >=9 ) ) {
3187 
3188  strncpy(result, line+strIndx,9) ;
3189  result[9] = '\0';
3190 
3191  Real localCharge = atof(result);
3192 
3193  // If we are reading charges from QM atoms, store them.
3194  if (atmIndx < msg->numQMAtoms ) {
3195  atmP[atmIndx].charge = localCharge;
3196  resForce[atmIndx].charge = localCharge;
3197  }
3198 
3199  // If we are reading charges from Dummy atoms,
3200  // place them on the appropriate QM atoms.
3201  if ( msg->numQMAtoms <= atmIndx &&
3202  atmIndx < msg->numAllAtoms ) {
3203  // We keep the calculated charge in the dummy atom
3204  // so that we can calculate electrostatic forces later on.
3205  atmP[atmIndx].charge = localCharge;
3206 
3207  // The dummy atom points to the QM atom to which it is bound.
3208  int qmInd = atmP[atmIndx].bountToIndx ;
3209  resForce[qmInd].charge += localCharge;
3210  }
3211 
3212  strIndx += 9;
3213  atmIndx++;
3214 
3215  // If we found all charges for QM and dummy atoms, break the loop
3216  // and stop reading the line.
3217  if (atmIndx == msg->numAllAtoms) {
3218  chargeFields = false;
3219  break;
3220  }
3221 
3222  }
3223 
3224  }
3225 
3226  int gradLength ; // Change for variable length MOPAC output
3227  if (gradFields) {
3228  if (atmIndx == 0) {
3229  double buf[10];
3230  int numfirstline = sscanf(line,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
3231  &buf[0],&buf[1],&buf[2],&buf[3],&buf[4],&buf[5],
3232  &buf[6],&buf[7],&buf[8],&buf[9]);
3233  gradLength = strlen(line)/numfirstline;
3234  }
3235  while ((strIndx < (strlen(line)-gradLength)) && (strlen(line)-1 >= gradLength ) ) {
3236 
3237  strncpy(result, line+strIndx,gradLength) ;
3238  result[gradLength] = '\0';
3239 
3240  gradient[gradCount] = atof(result);
3241  if (gradCount == 2) {
3242 
3243  if (atmIndx < msg->numQMAtoms ) {
3244  // Gradients in MOPAC are written in kcal/mol/A.
3245  resForce[atmIndx].force.x = -1*gradient[0];
3246  resForce[atmIndx].force.y = -1*gradient[1];
3247  resForce[atmIndx].force.z = -1*gradient[2];
3248  }
3249 
3250  // If we are reading forces applied on Dummy atoms,
3251  // place them on the appropriate QM and MM atoms to conserve energy.
3252 
3253  // This implementation was based on the description in
3254  // DOI: 10.1002/jcc.20857 : Equations 30 to 32
3255  if ( msg->numQMAtoms <= atmIndx &&
3256  atmIndx < msg->numAllAtoms ) {
3257  // The dummy atom points to the QM atom to which it is bound.
3258  // The QM atom and the MM atom (in a QM-MM bond) point to each other.
3259  int qmInd = atmP[atmIndx].bountToIndx ;
3260  int mmInd = atmP[qmInd].bountToIndx ;
3261 
3262  Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
3263  Real mmqmDist = dir.length() ;
3264 
3265  Real linkDist = Vector(atmP[atmIndx].position -
3266  atmP[qmInd].position).length() ;
3267 
3268  Force mmForce(0), qmForce(0),
3269  linkForce(gradient[0], gradient[1], gradient[2]);
3270  linkForce *= -1;
3271 
3272  Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
3273  // Unit vectors
3274  Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
3275  Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
3276  Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
3277  Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
3278 
3279  qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv +
3280  (xDelta)*base) )*xuv;
3281 
3282  qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv +
3283  (yDelta)*base) )*yuv;
3284 
3285  qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv +
3286  (zDelta)*base) )*zuv;
3287 
3288 
3289  mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
3290  (xDelta)*base) )*xuv;
3291 
3292  mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
3293  (yDelta)*base) )*yuv;
3294 
3295  mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
3296  (zDelta)*base) )*zuv;
3297 
3298  resForce[qmInd].force += qmForce;
3299  resForce[msg->numQMAtoms + mmInd].force += mmForce;
3300  }
3301 
3302  gradCount = 0;
3303  atmIndx++;
3304  } else {
3305  gradCount++;
3306  }
3307 
3308  strIndx += gradLength;
3309 
3310  // If we found all gradients for QM atoms, break the loop
3311  // and keep the next gradient line from being read, if there
3312  // is one. Following gradients, if present, will refer to link
3313  // atoms, and who cares about those?.
3314  if (atmIndx == msg->numAllAtoms) {
3315  gradFields = false;
3316  break;
3317  }
3318  }
3319  }
3320 
3321  }
3322 
3323  delete [] line;
3324 
3325  fclose(outputFile);
3326 
3327  // In case charges are not to be read form the QM software output,
3328  // we load the origianl atom charges.
3329  if (msg->qmAtmChrgMode == QMCHRGNONE) {
3330  int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
3331  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
3332  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
3333 
3334  atmIndx = 0 ;
3335  for (; atmIndx < msg->numQMAtoms; atmIndx++) {
3336 
3337  // Loops over all QM atoms (in all QM groups) comparing their global indices
3338  for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
3339 
3340  if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
3341 
3342  atmP[atmIndx].charge = qmAtmChrg[qmIter];
3343  resForce[atmIndx].charge = qmAtmChrg[qmIter];
3344 
3345  break;
3346  }
3347 
3348  }
3349 
3350  }
3351  for (size_t j=msg->numQMAtoms; j<msg->numAllAtoms; ++j ) {
3352  atmP[j].charge = 0;
3353  }
3354  }
3355 
3356  // remove force file
3357 // DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
3358 // iret = remove(outputFileName);
3359 // if ( iret ) { NAMD_err(0); }
3360 
3361  if (! (chargesRead && gradsRead) ) {
3362  NAMD_die("Error reading QM forces file. Not all data could be read!");
3363  }
3364 
3365  DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
3366 
3367  atmP = msg->data ;
3368  pcP = msg->data + msg->numAllAtoms ;
3369 
3370  // The initial point charge index for the *force* message is the number of QM
3371  // atoms, since the dummy atoms have no representation in NAMD
3372  int pcIndx = msg->numQMAtoms;
3373 
3374  // ---- WARNING ----
3375  // We need to apply the force felt by each QM atom due to the classical
3376  // point charges sent to MOPAC.
3377  // MOPAC takes the MM electrostatic potential into cosideration ONLY
3378  // when performing the SCF calculation and when returning the energy
3379  // of the system, but it does not apply the potential to the final
3380  // gradient calculation, therefore, we calculate the resulting force
3381  // here.
3382  // Initialize force vector for electrostatic contribution
3383  std::vector<Force> qmElecForce ;
3384  for (size_t j=0; j<msg->numAllAtoms; ++j )
3385  qmElecForce.push_back( Force(0) );
3386 
3387  // We loop over point charges and distribute the forces applied over
3388  // virtual point charges to the MM1 and MM2 atoms (if any virtual PCs exists).
3389  for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
3390 
3391  // No force was applied to the QM region due to this charge, so it
3392  // does not receive any back from the QM region. It must be an MM1 atom
3393  // from a QM-MM bond.
3394  if (pcP[i].type == QMPCTYPE_IGNORE)
3395  continue;
3396 
3397  Force totalForce(0);
3398 
3399  BigReal pntCharge = pcP[i].charge;
3400 
3401  Position posMM = pcP[i].position ;
3402 
3403  for (size_t j=0; j<msg->numAllAtoms; ++j ) {
3404 
3405  BigReal qmCharge = atmP[j].charge ;
3406 
3407  BigReal force = pntCharge*qmCharge*constants ;
3408 
3409  Position rVec = posMM - atmP[j].position ;
3410 
3411  force /= rVec.length2();
3412 
3413  // We accumulate the total force felt by a point charge
3414  // due to the QM charge distribution. This total force
3415  // will be applied to the point charge if it is a real one,
3416  // or will be distirbuted over MM1 and MM2 point charges, it
3417  // this is a virtual point charge.
3418  totalForce += force*rVec.unit();
3419 
3420  // Accumulate force felt by a QM atom due to point charges.
3421  qmElecForce[j] += -1*force*rVec.unit();
3422  }
3423 
3424  if (pcP[i].type == QMPCTYPE_CLASSICAL) {
3425  // Checking pcP was not a QM atom in another region
3426  // If so the interaction is accounted in the other region
3427  if (qmAtomGroup[pcP[i].id] == 0) {
3428  // If we already ignored MM1 charges, we take all other
3429  // non-virtual charges and apply forces directly to them.
3430  resForce[pcIndx].force += totalForce;
3431  }
3432  }
3433  else {
3434  // If we are handling virtual PC, we distribute the force over
3435  // MM1 and MM2.
3436 
3437  // Virtual PC are bound to MM2.
3438  int mm2Indx = pcP[i].bountToIndx;
3439  // MM2 charges are bound to MM1.
3440  int mm1Indx = pcP[mm2Indx].bountToIndx;
3441 
3442  Real Cq = pcP[i].dist;
3443 
3444  Force mm1Force = (1-Cq)*totalForce ;
3445  Force mm2Force = Cq*totalForce ;
3446 
3447  resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
3448  resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
3449  }
3450 
3451  }
3452 
3453  // We now apply the accumulated electrostatic force contribution
3454  // to QM atoms.
3455  for (size_t j=0; j<msg->numAllAtoms; ++j ) {
3456 
3457  if (j < msg->numQMAtoms ) {
3458  resForce[j].force += qmElecForce[j];
3459 
3460  } else {
3461  // In case the QM atom is a Link atom, we redistribute
3462  // its force as bove.
3463 
3464  // The dummy atom points to the QM atom to which it is bound.
3465  // The QM atom and the MM atom (in a QM-MM bond) point to each other.
3466  int qmInd = atmP[j].bountToIndx ;
3467  int mmInd = atmP[qmInd].bountToIndx ;
3468 
3469  Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
3470  Real mmqmDist = dir.length() ;
3471 
3472  Real linkDist = Vector(atmP[j].position -
3473  atmP[qmInd].position).length() ;
3474 
3475  Force linkForce = qmElecForce[j];
3476 
3477  Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
3478  // Unit vectors
3479  Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
3480  Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
3481  Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
3482  Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
3483 
3484  Force qmForce = (linkForce*((1 - linkDist/mmqmDist)*xuv +
3485  (xDelta)*base) )*xuv;
3486 
3487  qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv +
3488  (yDelta)*base) )*yuv;
3489 
3490  qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv +
3491  (zDelta)*base) )*zuv;
3492 
3493 
3494  Force mmForce = (linkForce*((linkDist/mmqmDist)*xuv -
3495  (xDelta)*base) )*xuv;
3496 
3497  mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
3498  (yDelta)*base) )*yuv;
3499 
3500  mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
3501  (zDelta)*base) )*zuv;
3502 
3503  resForce[qmInd].force += qmForce;
3504  resForce[msg->numQMAtoms + mmInd].force += mmForce;
3505  }
3506  }
3507 
3508  // Adjusts forces from PME, canceling contributions from the QM and
3509  // direct Coulomb forces calculated here.
3510  if (msg->PMEOn) {
3511 
3512  DebugM(1,"Correcting forces and energy for PME.\n");
3513 
3514  ewaldcof = msg->PMEEwaldCoefficient;
3515  BigReal TwoBySqrtPi = 1.12837916709551;
3516  pi_ewaldcof = TwoBySqrtPi * ewaldcof;
3517 
3518  for (size_t i=0; i < msg->numQMAtoms; i++) {
3519 
3520  BigReal p_i_charge = atmP[i].charge ;
3521  Position pos_i = atmP[i].position ;
3522 
3523  const BigReal kq_i = p_i_charge * constants;
3524 
3525  for (size_t j=i+1; j < msg->numQMAtoms; j++) {
3526 
3527  BigReal p_j_charge = atmP[j].charge ;
3528 
3529  Position pos_j = atmP[j].position ;
3530 
3531  BigReal r = Vector(pos_i - pos_j).length() ;
3532 
3533  BigReal tmp_a = r * ewaldcof;
3534  BigReal tmp_b = erfc(tmp_a);
3535  BigReal corr_energy = tmp_b;
3536  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
3537 
3538 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
3539  BigReal recip_energy = (1-tmp_b)/r;
3540 
3541 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
3542  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
3543 
3544  // Final force and energy correction for this pair of atoms.
3545  BigReal energy = kq_i * p_j_charge * recip_energy ;
3546 
3547  Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
3548 
3549  // The force is *subtracted* from the total force acting on
3550  // both atoms. The sign on fixForce corrects the orientation
3551  // of the subtracted force.
3552 // DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
3553 // << std::endl);
3554 // DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
3555 // << std::endl);
3556 // DebugM(4,"Force correction: " << fixForce << std::endl);
3557 // DebugM(4,"Energy correction: " << energy << "\n");
3558  resForce[i].force -= fixForce ;
3559  resForce[j].force -= -1*fixForce ;
3560 
3561  // The energy is *subtracted* from the total energy calculated here.
3562  resMsg->energyCorr -= energy;
3563  }
3564  }
3565 
3566 // DebugM(4,"Corrected energy QM-QM interactions: " << resMsg->energyCorr << "\n");
3567 
3568  pcIndx = msg->numQMAtoms;
3569  // We only loop over point charges from real atoms, ignoring the ones
3570  // created to handle QM-MM bonds.
3571  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
3572 
3573  BigReal p_i_charge = pcPpme[i].charge;
3574  Position pos_i = pcPpme[i].position ;
3575 
3576  const BigReal kq_i = p_i_charge * constants;
3577 
3578  Force fixForce = 0;
3579 
3580  for (size_t j=0; j<msg->numQMAtoms; ++j ) {
3581 
3582  BigReal p_j_charge = atmP[j].charge ;
3583 
3584  Position pos_j = atmP[j].position ;
3585 
3586  BigReal r = Vector(pos_i - pos_j).length() ;
3587 
3588  BigReal tmp_a = r * ewaldcof;
3589  BigReal tmp_b = erfc(tmp_a);
3590  BigReal corr_energy = tmp_b;
3591  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
3592 
3593 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
3594  BigReal recip_energy = (1-tmp_b)/r;
3595 
3596 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
3597  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
3598 
3599  // Final force and energy correction for this pair of atoms.
3600  BigReal energy = kq_i * p_j_charge * recip_energy ;
3601 
3602  fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
3603 
3604  // The energy is *subtracted* from the total energy calculated here.
3605  resMsg->energyCorr -= energy;
3606  }
3607 
3608  // The force is *subtracted* from the total force acting on
3609  // the point charge..
3610 // DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
3611 // << std::endl);
3612 // DebugM(4,"Force correction: " << fixForce << std::endl);
3613  resForce[pcIndx].force -= kq_i*fixForce ;
3614  }
3615 
3616  }
3617 
3618  // Calculates the virial contribution form the forces on QM atoms and
3619  // point charges calculated here.
3620  for (size_t i=0; i < msg->numQMAtoms; i++) {
3621  // virial used by NAMD is -'ve of normal convention, so reverse it!
3622  // virial[i][j] in file should be sum of -1 * f_i * r_j
3623  for ( int k=0; k<3; ++k )
3624  for ( int l=0; l<3; ++l )
3625  resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
3626  }
3627 
3628  pcIndx = msg->numQMAtoms; // Index in the real PC force array.
3629  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
3630  // virial used by NAMD is -'ve of normal convention, so reverse it!
3631  // virial[i][j] in file should be sum of -1 * f_i * r_j
3632  for ( int k=0; k<3; ++k )
3633  for ( int l=0; l<3; ++l )
3634  resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
3635  }
3636 
3637 
3638 
3639  // Send message to rank zero with results.
3640  QMProxy[0].recvQMRes(resMsg);
3641 
3642  if (msg->switching && msg->PMEOn)
3643  delete [] pcPpme;
3644 
3645  delete msg;
3646  return ;
3647 }
Node::Object
static Node * Object()
Definition: Node.h:86
QMAtomData::bountToIndx
int bountToIndx
Definition: ComputeQM.C:299
QMGrpCalcMsg::PMEEwaldCoefficient
BigReal PMEEwaldCoefficient
Definition: ComputeQM.C:336
QMGrpResMsg::force
QMForce * force
Definition: ComputeQM.C:181
QMGrpCalcMsg::baseDir
char baseDir[256]
Definition: ComputeQM.C:338
QMGrpCalcMsg::execPath
char execPath[256]
Definition: ComputeQM.C:338
NAMD_err
void NAMD_err(const char *err_msg)
Definition: common.C:106
QMGrpCalcMsg::data
QMAtomData * data
Definition: ComputeQM.C:339
Molecule::get_qmAtmIndx
const int * get_qmAtmIndx()
Definition: Molecule.h:824
QMGrpCalcMsg::constants
BigReal constants
Definition: ComputeQM.C:328
QMAtomData::charge
float charge
Definition: ComputeQM.C:297
Vector
Definition: Vector.h:64
QMGrpCalcMsg::prepProcOn
bool prepProcOn
Definition: ComputeQM.C:330
Molecule::get_numQMAtoms
int get_numQMAtoms()
Definition: Molecule.h:826
Real
float Real
Definition: common.h:109
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
QMGrpCalcMsg::numAllAtoms
int numAllAtoms
Definition: ComputeQM.C:324
QMPCTYPE_CLASSICAL
#define QMPCTYPE_CLASSICAL
Definition: ComputeQM.C:198
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
Vector::z
BigReal z
Definition: Vector.h:66
QMAtomData::dist
Real dist
Definition: ComputeQM.C:302
QMAtomData::position
Position position
Definition: ComputeQM.C:296
iout
#define iout
Definition: InfoStream.h:51
Vector::length
BigReal length(void) const
Definition: Vector.h:169
QMGrpCalcMsg::switching
bool switching
Definition: ComputeQM.C:332
Molecule::get_qmAtomGroup
const Real * get_qmAtomGroup() const
Definition: Molecule.h:820
QMGrpResMsg::numForces
int numForces
Definition: ComputeQM.C:180
QMGrpResMsg::energyOrig
BigReal energyOrig
Definition: ComputeQM.C:177
QMGrpCalcMsg::numQMAtoms
int numQMAtoms
Definition: ComputeQM.C:323
QMGrpCalcMsg::qmAtmChrgMode
int qmAtmChrgMode
Definition: ComputeQM.C:337
QMPCTYPE_IGNORE
#define QMPCTYPE_IGNORE
Definition: ComputeQM.C:197
z
gridSize z
Definition: ComputePmeCUDAKernel.cu:431
Bool
int Bool
Definition: common.h:133
QMGrpCalcMsg::numRealPntChrgs
int numRealPntChrgs
Definition: ComputeQM.C:325
Molecule::get_qmAtmChrg
Real * get_qmAtmChrg()
Definition: Molecule.h:823
QMGrpCalcMsg::prepProc
char prepProc[256]
Definition: ComputeQM.C:338
Vector::x
BigReal x
Definition: Vector.h:66
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:85
QMAtomData::element
char element[3]
Definition: ComputeQM.C:301
Vector::length2
BigReal length2(void) const
Definition: Vector.h:173
QMGrpCalcMsg::secProc
char secProc[256]
Definition: ComputeQM.C:338
Vector::y
BigReal y
Definition: Vector.h:66
charge
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
Definition: ComputeNonbondedBase2.h:179
y
gridSize y
Definition: ComputePmeCUDAKernel.cu:430
iERROR
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:83
QMGrpCalcMsg::numAllPntChrgs
int numAllPntChrgs
Definition: ComputeQM.C:326
x
gridSize x
Definition: ComputePmeCUDAKernel.cu:429
QMGrpResMsg::virial
BigReal virial[3][3]
Definition: ComputeQM.C:179
QMGrpCalcMsg::configLines
char * configLines
Definition: ComputeQM.C:340
Node::molecule
Molecule * molecule
Definition: Node.h:176
QMGrpCalcMsg::secProcOn
bool secProcOn
Definition: ComputeQM.C:329
Vector::unit
Vector unit(void) const
Definition: Vector.h:182
BigReal
double BigReal
Definition: common.h:114
QMCHRGNONE
#define QMCHRGNONE
Definition: Molecule.h:133
QMGrpResMsg::energyCorr
BigReal energyCorr
Definition: ComputeQM.C:178
for
for(int i=0;i< n1;++i)
Definition: ComputeFullDirectBase.h:34
void ComputeQMMgr::calcORCA ( QMGrpCalcMsg *  msg)

Definition at line 3649 of file ComputeQM.C.

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

3650 {
3651  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
3652  FILE *inputFile,*outputFile,*chrgFile;
3653  int iret;
3654 
3655  const size_t lineLen = 256;
3656  char *line = new char[lineLen];
3657 
3658  std::string qmCommand, inputFileName, outputFileName, pntChrgFileName, pcGradFilename;
3659  std::string tmpRedirectFileName, pcGradFileName;
3660 
3661  // For coulomb calculation
3662  BigReal constants = msg->constants;
3663 
3664  double gradient[3];
3665 
3666  DebugM(4,"Running ORCA on PE " << CkMyPe() << std::endl);
3667 
3668  QMAtomData *pcP = msg->data + msg->numAllAtoms ;
3669 
3670  // We create a pointer for PME correction, which may point to
3671  // a copy of the original point charge array, with unchanged charges,
3672  // or points to the original array in case no switching or charge
3673  // scheme is being used.
3674  QMAtomData *pcPpme = NULL;
3675  if (msg->switching) {
3676 
3677  if (msg->PMEOn)
3678  pcPpme = new QMAtomData[msg->numRealPntChrgs];
3679 
3680  pntChrgSwitching(msg, pcPpme) ;
3681  } else {
3682  pcPpme = pcP;
3683  }
3684 
3685  int retVal = 0;
3686  struct stat info;
3687 
3688  // For each QM group, create a subdirectory where files will be palced.
3689  std::string baseDir(msg->baseDir);
3690  std::ostringstream itosConv ;
3691  if ( CmiNumPartitions() > 1 ) {
3692  baseDir.append("/") ;
3693  itosConv << CmiMyPartition() ;
3694  baseDir += itosConv.str() ;
3695  itosConv.str("");
3696  itosConv.clear() ;
3697 
3698  if (stat(msg->baseDir, &info) != 0 ) {
3699  CkPrintf( "Node %d cannot access directory %s\n",
3700  CkMyPe(), baseDir.c_str() );
3701  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
3702  }
3703  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
3704  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
3705  retVal = mkdir(baseDir.c_str(), S_IRWXU);
3706  }
3707 
3708  }
3709  baseDir.append("/") ;
3710  itosConv << msg->grpIndx ;
3711  baseDir += itosConv.str() ;
3712 
3713  if (stat(msg->baseDir, &info) != 0 ) {
3714  CkPrintf( "Node %d cannot access directory %s\n",
3715  CkMyPe(), baseDir.c_str() );
3716  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
3717  }
3718  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
3719  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
3720  int retVal = mkdir(baseDir.c_str(), S_IRWXU);
3721  }
3722 
3723  #ifdef DEBUG_QM
3724  Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
3725  #endif
3726 
3727  inputFileName.clear();
3728  inputFileName.append(baseDir.c_str()) ;
3729  inputFileName.append("/qmmm_") ;
3730  inputFileName += itosConv.str() ;
3731  inputFileName.append(".input") ;
3732 
3733  // Opens file for coordinate and parameter input
3734  inputFile = fopen(inputFileName.c_str(),"w");
3735  if ( ! inputFile ) {
3736  iout << iERROR << "Could not open input file for writing: "
3737  << inputFileName << "\n" << endi ;
3738  NAMD_err("Error writing QM input file.");
3739  }
3740 
3741  // Builds the command that will be executed
3742  qmCommand.clear();
3743  qmCommand.append("cd ");
3744  qmCommand.append(baseDir);
3745  qmCommand.append(" ; ");
3746  qmCommand.append(msg->execPath) ;
3747  qmCommand.append(" ") ;
3748  qmCommand.append(inputFileName) ;
3749 
3750  // Build the redirect file name we need for the screen output.
3751  // That's the only place where we can get partial charges for QM atoms.
3752  tmpRedirectFileName = inputFileName ;
3753  tmpRedirectFileName.append(".TmpOut") ;
3754 
3755  qmCommand.append(" > ") ;
3756  qmCommand.append(tmpRedirectFileName) ;
3757 
3758  // Builds the file name where orca will place the gradient
3759  // This will be relative to the input file
3760  outputFileName = inputFileName ;
3761  outputFileName.append(".engrad") ;
3762 
3763  // Builds the file name where orca will place gradients acting on
3764  // point charges
3765  pcGradFilename = inputFileName ;
3766  pcGradFilename.append(".pcgrad") ;
3767 
3768  if (msg->numAllPntChrgs) {
3769  // Builds the file name where we will write the point charges.
3770  pntChrgFileName = inputFileName ;
3771  pntChrgFileName.append(".pntchrg") ;
3772 
3773  pcGradFileName = inputFileName;
3774  pcGradFileName.append(".pcgrad");
3775 
3776  chrgFile = fopen(pntChrgFileName.c_str(),"w");
3777  if ( ! chrgFile ) {
3778  iout << iERROR << "Could not open charge file for writing: "
3779  << pntChrgFileName << "\n" << endi ;
3780  NAMD_err("Error writing point charge file.");
3781  }
3782 
3783  int numPntChrgs = 0;
3784  for (int i=0; i<msg->numAllPntChrgs; i++ ) {
3785  if (pcP[i].type != QMPCTYPE_IGNORE)
3786  numPntChrgs++;
3787  }
3788 
3789  iret = fprintf(chrgFile,"%d\n", numPntChrgs);
3790  if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
3791  }
3792 
3793  // writes configuration lines to the input file.
3794  std::stringstream ss(msg->configLines) ;
3795  std::string confLineStr;
3796  while (std::getline(ss, confLineStr) ) {
3797  confLineStr.append("\n");
3798  iret = fprintf(inputFile,confLineStr.c_str());
3799  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3800  }
3801 
3802  if (msg->numAllPntChrgs) {
3803  iret = fprintf(inputFile,"%%pointcharges \"%s\"\n", pntChrgFileName.c_str());
3804  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3805  }
3806 
3807  iret = fprintf(inputFile,"\n\n%%coords\n CTyp xyz\n");
3808  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3809 
3810  iret = fprintf(inputFile," Charge %f\n",msg->charge);
3811  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3812 
3813  iret = fprintf(inputFile," Mult %f\n",msg->multiplicity);
3814  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3815 
3816  iret = fprintf(inputFile," Units Angs\n coords\n\n");
3817  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3818 
3819  DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " <<
3820  inputFileName.c_str() << " and " << msg->numAllPntChrgs <<
3821  " point charges in file " << pntChrgFileName.c_str() << "\n");
3822 
3823  // write QM and dummy atom coordinates to input file.
3824  QMAtomData *atmP = msg->data ;
3825  for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
3826 
3827  double x = atmP->position.x;
3828  double y = atmP->position.y;
3829  double z = atmP->position.z;
3830 
3831  iret = fprintf(inputFile," %s %f %f %f\n",
3832  atmP->element,x,y,z);
3833  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3834 
3835  }
3836 
3837  iret = fprintf(inputFile," end\nend\n");
3838  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
3839 
3840  if (msg->numAllPntChrgs) {
3841  // Write point charges to file.
3842  pcP = msg->data + msg->numAllAtoms ;
3843  for ( size_t j=0; j < msg->numAllPntChrgs; j++, ++pcP) {
3844 
3845  if (pcP->type == QMPCTYPE_IGNORE)
3846  continue;
3847 
3848  double charge = pcP->charge;
3849 
3850  double x = pcP->position.x;
3851  double y = pcP->position.y;
3852  double z = pcP->position.z;
3853 
3854  iret = fprintf(chrgFile,"%f %f %f %f\n",
3855  charge,x,y,z);
3856  if ( iret < 0 ) { NAMD_err("Error writing point charge file."); }
3857  }
3858 
3859  fclose(chrgFile);
3860  }
3861 
3862  DebugM(4,"Closing input and charge file\n");
3863  fclose(inputFile);
3864 
3865  if (msg->prepProcOn) {
3866 
3867  std::string prepProc(msg->prepProc) ;
3868  prepProc.append(" ") ;
3869  prepProc.append(inputFileName) ;
3870  iret = system(prepProc.c_str());
3871  if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
3872  if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
3873  }
3874 
3875  // runs QM command
3876  DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
3877  iret = system(qmCommand.c_str());
3878 
3879  if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
3880  if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
3881 
3882  if (msg->secProcOn) {
3883 
3884  std::string secProc(msg->secProc) ;
3885  secProc.append(" ") ;
3886  secProc.append(inputFileName) ;
3887  itosConv.str("");
3888  itosConv.clear() ;
3889  itosConv << msg->timestep ;
3890  secProc.append(" ") ;
3891  secProc += itosConv.str() ;
3892 
3893  iret = system(secProc.c_str());
3894  if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
3895  if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
3896  }
3897 
3898  // remove coordinate file
3899 // iret = remove(inputFileName);
3900 // if ( iret ) { NAMD_err(0); }
3901 
3902  // remove coordinate file
3903 // iret = remove(pntChrgFileName);
3904 // if ( iret ) { NAMD_err(0); }
3905 
3906  // opens output file
3907  DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
3908  outputFile = fopen(outputFileName.c_str(),"r");
3909  if ( ! outputFile ) {
3910  iout << iERROR << "Could not find QM output file!\n" << endi;
3911  NAMD_err(0);
3912  }
3913 
3914  // Resets the pointers.
3915  atmP = msg->data ;
3916  pcP = msg->data + msg->numAllAtoms ;
3917 
3918  // Allocates the return message.
3919  QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
3920  resMsg->grpIndx = msg->grpIndx;
3921  resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
3922  resMsg->energyOrig = 0;
3923  resMsg->energyCorr = 0;
3924  for ( int k=0; k<3; ++k )
3925  for ( int l=0; l<3; ++l )
3926  resMsg->virial[k][l] = 0;
3927  QMForce *resForce = resMsg->force ;
3928 
3929  // We split the data into two pointers so we can "skip" the dummy atoms
3930  // which have no representation in NAMD.
3931  for (int i=0; i<resMsg->numForces; i++) {
3932  resForce[i].force = 0;
3933  resForce[i].charge = 0 ;
3934  if (i < msg->numQMAtoms) {
3935  // We use the replace field to indicate QM atoms,
3936  // which will have charge information.
3937  resForce[i].replace = 1 ;
3938  resForce[i].id = atmP->id;
3939  atmP++;
3940  }
3941  else {
3942  // We use the replace field to indicate QM atoms,
3943  // which will have charge information.
3944  resForce[i].replace = 0 ;
3945  resForce[i].id = pcP->id;
3946  pcP++;
3947  }
3948  }
3949 
3950  // Resets the pointers.
3951  atmP = msg->data ;
3952  pcP = msg->data + msg->numAllAtoms ;
3953 
3954  size_t atmIndx = 0, gradCount = 0;
3955  int numAtomsInFile = 0;
3956 
3957  // Reads the data form the output file created by the QM software.
3958  // Gradients over the QM atoms, and Charges for QM atoms will be read.
3959 
3960  // skip lines before number of atoms
3961  for (int i = 0; i < 3; i++) {
3962  fgets(line, lineLen, outputFile);
3963  }
3964 
3965  iret = fscanf(outputFile,"%d\n", &numAtomsInFile);
3966  if ( iret != 1 ) {
3967  NAMD_die("Error reading QM forces file.");
3968  }
3969 
3970  #ifdef DEBUG_QM
3971  if(numAtomsInFile != msg->numAllAtoms) {
3972  NAMD_die("Error reading QM forces file. Number of atoms in QM output\
3973  does not match the expected.");
3974  }
3975  #endif
3976 
3977  // skip lines before energy
3978  for (int i = 0; i < 3; i++) {
3979  fgets(line, lineLen, outputFile);
3980  }
3981 
3982  iret = fscanf(outputFile,"%lf\n", &resMsg->energyOrig);
3983  if ( iret != 1 ) {
3984  NAMD_die("Error reading QM forces file.");
3985  }
3986 // iout << "Energy in step (Hartree): " << resMsg->energyOrig << "\n" << endi;
3987  // All energies are given in Eh (Hartree)
3988  // NAMD needs energies in kcal/mol
3989  // The conversion factor is 627.509469
3990  resMsg->energyOrig *= 627.509469;
3991 // iout << "Energy in step (Kcal/mol): " << resMsg->energyOrig << "\n" << endi;
3992 
3993  resMsg->energyCorr = resMsg->energyOrig;
3994 
3995  // skip lines before gradient
3996  for (int i = 0; i < 3; i++) {
3997  fgets(line, lineLen, outputFile) ;
3998  }
3999 
4000  // Break the loop when we find all gradients for QM atoms,
4001  // and keep the next gradient lines from being read, if there
4002  // are more. Following gradients, if present, will refer to link
4003  // atoms.
4004  atmIndx = 0;
4005  gradCount = 0;
4006  for ( size_t i=0; i<msg->numAllAtoms*3; ++i ) {
4007 
4008  iret = fscanf(outputFile,"%lf\n", &gradient[gradCount]);
4009  if ( iret != 1 ) { NAMD_die("Error reading QM forces file."); }
4010 
4011  if (gradCount == 2){
4012 
4013  // All gradients are given in Eh/a0 (Hartree over Bohr radius)
4014  // NAMD needs forces in kcal/mol/angstrons
4015  // The conversion factor is 627.509469/0.529177 = 1185.82151
4016  if (atmIndx < msg->numQMAtoms ) {
4017  resForce[atmIndx].force.x = -1*gradient[0]*1185.82151;
4018  resForce[atmIndx].force.y = -1*gradient[1]*1185.82151;
4019  resForce[atmIndx].force.z = -1*gradient[2]*1185.82151;
4020  }
4021 
4022  // If we are reading forces applied on Dummy atoms,
4023  // place them on the appropriate QM and MM atoms to conserve energy.
4024 
4025  // This implementation was based on the description in
4026  // DOI: 10.1002/jcc.20857 : Equations 30 to 32
4027  if ( msg->numQMAtoms <= atmIndx &&
4028  atmIndx < msg->numAllAtoms ) {
4029  // The dummy atom points to the QM atom to which it is bound.
4030  // The QM atom and the MM atom (in a QM-MM bond) point to each other.
4031  int qmInd = atmP[atmIndx].bountToIndx ;
4032  int mmInd = atmP[qmInd].bountToIndx ;
4033 
4034  Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
4035  Real mmqmDist = dir.length() ;
4036 
4037  Real linkDist = Vector(atmP[atmIndx].position - atmP[qmInd].position).length() ;
4038 
4039  Force mmForce(0), qmForce(0), linkForce(gradient[0], gradient[1], gradient[2]);
4040  linkForce *= -1*1185.82151;
4041 
4042  Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
4043  // Unit vectors
4044  Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
4045  Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
4046  Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
4047  Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
4048 
4049  qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv +
4050  (xDelta)*base) )*xuv;
4051 
4052  qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv +
4053  (yDelta)*base) )*yuv;
4054 
4055  qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv +
4056  (zDelta)*base) )*zuv;
4057 
4058 
4059  mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
4060  (xDelta)*base) )*xuv;
4061 
4062  mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
4063  (yDelta)*base) )*yuv;
4064 
4065  mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
4066  (zDelta)*base) )*zuv;
4067 
4068  resForce[qmInd].force += qmForce;
4069  resForce[msg->numQMAtoms + mmInd].force += mmForce;
4070  }
4071 
4072  gradCount = 0;
4073  atmIndx++ ;
4074  } else
4075  gradCount++ ;
4076 
4077  }
4078 
4079  fclose(outputFile);
4080 
4081  // In case charges are not to be read form the QM software output,
4082  // we load the origianl atom charges.
4083  if (msg->qmAtmChrgMode == QMCHRGNONE) {
4084  int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
4085  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
4086  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
4087 
4088  atmIndx = 0 ;
4089  for (; atmIndx < msg->numQMAtoms; atmIndx++) {
4090 
4091  // Loops over all QM atoms (in all QM groups) comparing their global indices
4092  for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
4093 
4094  if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
4095 
4096  atmP[atmIndx].charge = qmAtmChrg[qmIter];
4097  resForce[atmIndx].charge = qmAtmChrg[qmIter];
4098 
4099  break;
4100  }
4101 
4102  }
4103 
4104  }
4105  }
4106  else {
4107  // opens redirected output file
4108  outputFile = fopen(tmpRedirectFileName.c_str(),"r");
4109  if ( ! outputFile ) {
4110  iout << iERROR << "Could not find Redirect output file:"
4111  << tmpRedirectFileName << std::endl << endi;
4112  NAMD_err(0);
4113  }
4114 
4115  DebugM(4,"Opened tmeporary output for charge reading: " << tmpRedirectFileName.c_str() << "\n");
4116 
4117  int lineState = 0;
4118  char result[20] ;
4119  while ( fgets(line, lineLen, outputFile) != NULL){
4120 
4121  // We loop over the lines of the output file untill we find
4122  // the line that initiates the charges lines. Then
4123  // we read all lines and expect to get one charge
4124  // per line, untill we find a line that sums all charges.
4125 
4126  switch (msg->qmAtmChrgMode) {
4127 
4128  case QMCHRGMULLIKEN:
4129  {
4130 
4131  if ( strstr(line,"MULLIKEN ATOMIC CHARGES") != NULL ) {
4132  lineState = 1;
4133  atmIndx = 0;
4134 
4135  // Now we expect the following line to have a series of dashes
4136  // and the folowing lines to have atom charges (among other info)
4137  continue;
4138  }
4139 
4140  if ( strstr(line,"Sum of atomic charges") != NULL ) {
4141 
4142  // Now that all charges have been read, grabs the total charge
4143  // just for fun... and checking purposes.
4144  #ifdef DEBUG_QM
4145  strncpy(result, line + 31,12) ;
4146  result[12] = '\0';
4147 
4148  DebugM(4,"Total charge of QM region calculated by ORCA is: "
4149  << atof(result) << std::endl )
4150  #endif
4151 
4152  // Nothing more to read from the output file
4153  break;
4154  }
4155 
4156  // Line state 1 means we have to skip ONLY one line.
4157  if (lineState == 1) {
4158  lineState = 2;
4159  continue;
4160  }
4161 
4162  // Line state 2 means we have to read the line and grab the info.
4163  if (lineState == 2) {
4164 
4165  strncpy(result, line+8,12) ;
4166  result[12] = '\0';
4167 
4168  Real localCharge = atof(result);
4169 
4170  // If we are reading charges from QM atoms, store them.
4171  if (atmIndx < msg->numQMAtoms ) {
4172  atmP[atmIndx].charge = localCharge;
4173  resForce[atmIndx].charge = localCharge;
4174  }
4175 
4176  // If we are reading charges from Dummy atoms,
4177  // place the on the appropriate QM atom.
4178  if ( msg->numQMAtoms <= atmIndx &&
4179  atmIndx < msg->numAllAtoms ) {
4180  int qmInd = atmP[atmIndx].bountToIndx ;
4181  atmP[qmInd].charge += localCharge;
4182  resForce[qmInd].charge += localCharge;
4183  }
4184 
4185  atmIndx++ ;
4186 
4187  // If we found all charges for QM atoms, change the lineState
4188  // untill we reach the "total charge" line.
4189  if (atmIndx == msg->numAllAtoms ) {
4190  lineState = 0;
4191  }
4192 
4193  continue;
4194  }
4195 
4196  } break ;
4197 
4198  case QMCHRGCHELPG :
4199  {
4200 
4201  if ( strstr(line,"CHELPG Charges") != NULL ) {
4202  lineState = 1;
4203  atmIndx = 0;
4204 
4205  // Now we expect the following line to have a series of dashes
4206  // and the folowing lines to have atom charges (among other info)
4207  continue;
4208  }
4209 
4210  if ( strstr(line,"Total charge") != NULL ) {
4211 
4212  // Now that all charges have been read, grabs the total charge
4213  // just for fun... and checking purposes.
4214  #ifdef DEBUG_QM
4215  strncpy(result, line + 14,13) ;
4216  result[13] = '\0';
4217 
4218  DebugM(4,"Total charge of QM region calculated by ORCA is: "
4219  << atof(result) << std::endl )
4220  #endif
4221 
4222  // Nothing more to read from the output file
4223  break;
4224  }
4225 
4226  // Line state 1 means we have to skip ONLY one line.
4227  if (lineState == 1) {
4228  lineState = 2;
4229  continue;
4230  }
4231 
4232  // Line state 2 means we have to read the line and grab the info.
4233  if (lineState == 2) {
4234 
4235  strncpy(result, line+12,15) ;
4236  result[15] = '\0';
4237 
4238  Real localCharge = atof(result);
4239 
4240  // If we are reading charges from QM atoms, store them.
4241  if (atmIndx < msg->numQMAtoms ) {
4242  atmP[atmIndx].charge = localCharge;
4243  resForce[atmIndx].charge = localCharge;
4244  }
4245 
4246  // If we are reading charges from Dummy atoms,
4247  // place the on the appropriate QM atom.
4248  if ( msg->numQMAtoms <= atmIndx &&
4249  atmIndx < msg->numAllAtoms ) {
4250  int qmInd = atmP[atmIndx].bountToIndx ;
4251  atmP[qmInd].charge += localCharge;
4252  resForce[qmInd].charge += localCharge;
4253  }
4254 
4255  atmIndx++ ;
4256 
4257  // If we found all charges for QM atoms, we ignore the following line
4258  // untill we reach the "total charge" line.
4259  if (atmIndx == msg->numAllAtoms ) {
4260  lineState = 1;
4261  }
4262 
4263  continue;
4264  }
4265 
4266  } break;
4267  }
4268  }
4269 
4270  fclose(outputFile);
4271  }
4272 
4273  // remove force file
4274 // DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
4275 // iret = remove(outputFileName);
4276 // if ( iret ) { NAMD_err(0); }
4277 
4278 
4279  DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
4280 
4281  atmP = msg->data ;
4282  pcP = msg->data + msg->numAllAtoms ;
4283 
4284  // The initial point charge index for the force message is the number of QM
4285  // atoms, since the dummy atoms have no representation in NAMD
4286  int pcIndx = msg->numQMAtoms;
4287 
4288  // If there are no point charges, orca will not create a "pcgrad" file.
4289  if (msg->numAllPntChrgs > 0) {
4290 
4291  // opens redirected output file
4292  outputFile = fopen(pcGradFileName.c_str(),"r");
4293  if ( ! outputFile ) {
4294  iout << iERROR << "Could not find PC gradient output file:"
4295  << pcGradFileName << std::endl << endi;
4296  NAMD_err(0);
4297  }
4298 
4299  DebugM(4,"Opened pc output for gradient reading: " << pcGradFileName.c_str() << "\n");
4300 
4301  int pcgradNumPC = 0, readPC = 0;
4302  iret = fscanf(outputFile,"%d\n", &pcgradNumPC );
4303  if ( iret != 1 ) {
4304  NAMD_die("Error reading PC forces file.");
4305  }
4306  DebugM(4,"Found in pc gradient output: " << pcgradNumPC << " point charge grads.\n");
4307 
4308  // We loop over point charges, reading the total electrostatic force
4309  // applied on them by the QM region.
4310  // We redistribute the forces applied over virtual point
4311  // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
4312  for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
4313 
4314  Force totalForce(0);
4315 
4316  // No force was applied to the QM region due to this charge, since it
4317  // was skipped when writing down point charges to the QM software, so it
4318  // does not receive any back from the QM region. It must be an MM1 atom
4319  // from a QM-MM bond.
4320  if (pcP[i].type == QMPCTYPE_IGNORE)
4321  continue;
4322 
4323  fgets(line, lineLen, outputFile);
4324 
4325  iret = sscanf(line,"%lf %lf %lf\n", &totalForce[0], &totalForce[1], &totalForce[2] );
4326  if ( iret != 3 ) {
4327  NAMD_die("Error reading PC forces file.");
4328  }
4329  // All gradients are given in Eh/a0 (Hartree over Bohr radius)
4330  // NAMD needs forces in kcal/mol/angstrons
4331  // The conversion factor is 627.509469/0.529177 = 1185.82151
4332  totalForce *= -1185.82151;
4333 
4334  if (pcP[i].type == QMPCTYPE_CLASSICAL) {
4335  // Checking pcP was not a QM atom in another region
4336  // If so the interaction is accounted in the other region
4337  if (qmAtomGroup[pcP[i].id] == 0) {
4338  // If we already ignored MM1 charges, we take all other
4339  // non-virtual charges and apply forces directly to them.
4340  resForce[pcIndx].force += totalForce;
4341  }
4342  }
4343  else {
4344  // If we are handling virtual PC, we distribute the force over
4345  // MM1 and MM2.
4346 
4347  Force mm1Force(0), mm2Force(0);
4348 
4349  // Virtual PC are bound to MM2.
4350  int mm2Indx = pcP[i].bountToIndx;
4351  // MM2 charges are bound to MM1.
4352  int mm1Indx = pcP[mm2Indx].bountToIndx;
4353 
4354  Real Cq = pcP[i].dist;
4355 
4356  mm1Force = (1-Cq)*totalForce ;
4357  mm2Force = Cq*totalForce ;
4358 
4359  resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
4360  resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
4361  }
4362 
4363  }
4364 
4365  fclose(outputFile);
4366  }
4367 
4368  delete [] line;
4369 
4370  // Adjusts forces from PME, canceling contributions from the QM and
4371  // direct Coulomb forces calculated here.
4372  if (msg->PMEOn) {
4373 
4374  DebugM(1,"Correcting forces and energy for PME.\n");
4375 
4376  ewaldcof = msg->PMEEwaldCoefficient;
4377  BigReal TwoBySqrtPi = 1.12837916709551;
4378  pi_ewaldcof = TwoBySqrtPi * ewaldcof;
4379 
4380  for (size_t i=0; i < msg->numQMAtoms; i++) {
4381 
4382  BigReal p_i_charge = atmP[i].charge ;
4383  Position pos_i = atmP[i].position ;
4384 
4385  for (size_t j=i+1; j < msg->numQMAtoms; j++) {
4386 
4387  BigReal p_j_charge = atmP[j].charge ;
4388 
4389  Position pos_j = atmP[j].position ;
4390 
4391  BigReal r = Vector(pos_i - pos_j).length() ;
4392 
4393  BigReal tmp_a = r * ewaldcof;
4394  BigReal tmp_b = erfc(tmp_a);
4395  BigReal corr_energy = tmp_b;
4396  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
4397 
4398 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
4399  BigReal recip_energy = (1-tmp_b)/r;
4400 
4401 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
4402  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
4403 
4404  const BigReal kq_i = p_i_charge * constants;
4405 
4406  // Final force and energy correction for this pair of atoms.
4407  BigReal energy = kq_i * p_j_charge * recip_energy ;
4408 
4409  Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
4410 
4411  // The force is *subtracted* from the total force acting on
4412  // both atoms. The sign on fixForce corrects the orientation
4413  // of the subtracted force.
4414 // DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
4415 // << std::endl);
4416 // DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
4417 // << std::endl);
4418 // DebugM(4,"Force correction: " << fixForce << std::endl);
4419  resForce[i].force -= fixForce ;
4420  resForce[j].force -= -1*fixForce;
4421 
4422  // The energy is *subtracted* from the total energy calculated here.
4423  resMsg->energyCorr -= energy;
4424  }
4425  }
4426 
4427  pcIndx = msg->numQMAtoms;
4428  // We only loop over point charges from real atoms, ignoring the ones
4429  // created to handle QM-MM bonds.
4430  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
4431 
4432  BigReal p_i_charge = pcPpme[i].charge;
4433  Position pos_i = pcPpme[i].position ;
4434 
4435  const BigReal kq_i = p_i_charge * constants;
4436 
4437  Force fixForce = 0;
4438 
4439  for (size_t j=0; j<msg->numQMAtoms; ++j ) {
4440 
4441  BigReal p_j_charge = atmP[j].charge ;
4442 
4443  Position pos_j = atmP[j].position ;
4444 
4445  BigReal r = Vector(pos_i - pos_j).length() ;
4446 
4447  BigReal tmp_a = r * ewaldcof;
4448  BigReal tmp_b = erfc(tmp_a);
4449  BigReal corr_energy = tmp_b;
4450  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
4451 
4452 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
4453  BigReal recip_energy = (1-tmp_b)/r;
4454 
4455 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
4456  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
4457 
4458  // Final force and energy correction for this pair of atoms.
4459  BigReal energy = kq_i * p_j_charge * recip_energy ;
4460 
4461  fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
4462 
4463  // The energy is *subtracted* from the total energy calculated here.
4464  resMsg->energyCorr -= energy;
4465 
4466  }
4467 
4468  // The force is *subtracted* from the total force acting on
4469  // the point charge.
4470 // DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
4471 // << std::endl);
4472 // DebugM(4,"Force correction: " << fixForce << std::endl);
4473  resForce[pcIndx].force -= kq_i*fixForce ;
4474  }
4475 
4476  }
4477 
4478  DebugM(1,"Determining virial...\n");
4479 
4480  // Calculates the virial contribution form the forces on QM atoms and
4481  // point charges calculated here.
4482  for (size_t i=0; i < msg->numQMAtoms; i++) {
4483  // virial used by NAMD is -'ve of normal convention, so reverse it!
4484  // virial[i][j] in file should be sum of -1 * f_i * r_j
4485  for ( int k=0; k<3; ++k )
4486  for ( int l=0; l<3; ++l )
4487  resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
4488  }
4489 
4490  pcIndx = msg->numQMAtoms; // Index in the real PC force array.
4491  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
4492  // virial used by NAMD is -'ve of normal convention, so reverse it!
4493  // virial[i][j] in file should be sum of -1 * f_i * r_j
4494  for ( int k=0; k<3; ++k )
4495  for ( int l=0; l<3; ++l )
4496  resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
4497  }
4498 
4499  DebugM(1,"End of QM processing. Sending result message.\n");
4500 
4501  // Send message to rank zero with results.
4502  QMProxy[0].recvQMRes(resMsg);
4503 
4504  if (msg->switching && msg->PMEOn)
4505  delete [] pcPpme;
4506 
4507  delete msg;
4508  return ;
4509 }
Node::Object
static Node * Object()
Definition: Node.h:86
QMAtomData::bountToIndx
int bountToIndx
Definition: ComputeQM.C:299
QMGrpCalcMsg::PMEEwaldCoefficient
BigReal PMEEwaldCoefficient
Definition: ComputeQM.C:336
QMGrpResMsg::force
QMForce * force
Definition: ComputeQM.C:181
QMGrpCalcMsg::baseDir
char baseDir[256]
Definition: ComputeQM.C:338
QMGrpCalcMsg::execPath
char execPath[256]
Definition: ComputeQM.C:338
NAMD_err
void NAMD_err(const char *err_msg)
Definition: common.C:106
QMGrpCalcMsg::data
QMAtomData * data
Definition: ComputeQM.C:339
QMGrpCalcMsg::multiplicity
Real multiplicity
Definition: ComputeQM.C:327
Molecule::get_qmAtmIndx
const int * get_qmAtmIndx()
Definition: Molecule.h:824
QMGrpCalcMsg::constants
BigReal constants
Definition: ComputeQM.C:328
QMAtomData::charge
float charge
Definition: ComputeQM.C:297
Vector
Definition: Vector.h:64
QMGrpCalcMsg::prepProcOn
bool prepProcOn
Definition: ComputeQM.C:330
Molecule::get_numQMAtoms
int get_numQMAtoms()
Definition: Molecule.h:826
Real
float Real
Definition: common.h:109
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
QMGrpCalcMsg::numAllAtoms
int numAllAtoms
Definition: ComputeQM.C:324
QMPCTYPE_CLASSICAL
#define QMPCTYPE_CLASSICAL
Definition: ComputeQM.C:198
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
Vector::z
BigReal z
Definition: Vector.h:66
QMAtomData::dist
Real dist
Definition: ComputeQM.C:302
QMAtomData::position
Position position
Definition: ComputeQM.C:296
iout
#define iout
Definition: InfoStream.h:51
Vector::length
BigReal length(void) const
Definition: Vector.h:169
QMCHRGCHELPG
#define QMCHRGCHELPG
Definition: Molecule.h:135
QMGrpCalcMsg::switching
bool switching
Definition: ComputeQM.C:332
Molecule::get_qmAtomGroup
const Real * get_qmAtomGroup() const
Definition: Molecule.h:820
QMGrpResMsg::numForces
int numForces
Definition: ComputeQM.C:180
QMGrpResMsg::energyOrig
BigReal energyOrig
Definition: ComputeQM.C:177
QMGrpCalcMsg::numQMAtoms
int numQMAtoms
Definition: ComputeQM.C:323
QMGrpCalcMsg::qmAtmChrgMode
int qmAtmChrgMode
Definition: ComputeQM.C:337
QMCHRGMULLIKEN
#define QMCHRGMULLIKEN
Definition: Molecule.h:134
QMPCTYPE_IGNORE
#define QMPCTYPE_IGNORE
Definition: ComputeQM.C:197
z
gridSize z
Definition: ComputePmeCUDAKernel.cu:431
QMGrpCalcMsg::numRealPntChrgs
int numRealPntChrgs
Definition: ComputeQM.C:325
Molecule::get_qmAtmChrg
Real * get_qmAtmChrg()
Definition: Molecule.h:823
QMGrpCalcMsg::prepProc
char prepProc[256]
Definition: ComputeQM.C:338
Vector::x
BigReal x
Definition: Vector.h:66
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:85
QMGrpCalcMsg::secProc
char secProc[256]
Definition: ComputeQM.C:338
Vector::y
BigReal y
Definition: Vector.h:66
charge
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
Definition: ComputeNonbondedBase2.h:179
y
gridSize y
Definition: ComputePmeCUDAKernel.cu:430
iERROR
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:83
QMGrpCalcMsg::numAllPntChrgs
int numAllPntChrgs
Definition: ComputeQM.C:326
x
gridSize x
Definition: ComputePmeCUDAKernel.cu:429
QMGrpResMsg::virial
BigReal virial[3][3]
Definition: ComputeQM.C:179
QMGrpCalcMsg::configLines
char * configLines
Definition: ComputeQM.C:340
Node::molecule
Molecule * molecule
Definition: Node.h:176
QMGrpCalcMsg::secProcOn
bool secProcOn
Definition: ComputeQM.C:329
BigReal
double BigReal
Definition: common.h:114
QMCHRGNONE
#define QMCHRGNONE
Definition: Molecule.h:133
QMGrpResMsg::energyCorr
BigReal energyCorr
Definition: ComputeQM.C:178
for
for(int i=0;i< n1;++i)
Definition: ComputeFullDirectBase.h:34
void ComputeQMMgr::calcUSR ( QMGrpCalcMsg *  msg)

Definition at line 4511 of file ComputeQM.C.

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

4511  {
4512  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
4513  FILE *inputFile,*outputFile;
4514  int iret;
4515 
4516  std::string qmCommand, inputFileName, outputFileName;
4517 
4518  // For coulomb calculation
4519  BigReal constants = msg->constants;
4520 
4521  DebugM(4,"Running USER DEFINED SOFTWARE on PE " << CkMyPe() << std::endl);
4522 
4523  QMAtomData *pcP = msg->data + msg->numAllAtoms ;
4524 
4525  // We create a pointer for PME correction, which may point to
4526  // a copy of the original point charge array, with unchanged charges,
4527  // or points to the original array in case no switching or charge
4528  // scheme is being used.
4529  QMAtomData *pcPpme = NULL;
4530  if (msg->switching) {
4531 
4532  if (msg->PMEOn)
4533  pcPpme = new QMAtomData[msg->numRealPntChrgs];
4534 
4535  pntChrgSwitching(msg, pcPpme) ;
4536  } else {
4537  pcPpme = pcP;
4538  }
4539 
4540  int retVal = 0;
4541  struct stat info;
4542 
4543  // For each QM group, create a subdirectory where files will be palced.
4544  std::string baseDir(msg->baseDir);
4545  std::ostringstream itosConv ;
4546  if ( CmiNumPartitions() > 1 ) {
4547  baseDir.append("/") ;
4548  itosConv << CmiMyPartition() ;
4549  baseDir += itosConv.str() ;
4550  itosConv.str("");
4551  itosConv.clear() ;
4552 
4553  if (stat(msg->baseDir, &info) != 0 ) {
4554  CkPrintf( "Node %d cannot access directory %s\n",
4555  CkMyPe(), baseDir.c_str() );
4556  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
4557  }
4558  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
4559  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
4560  retVal = mkdir(baseDir.c_str(), S_IRWXU);
4561  }
4562 
4563  }
4564  baseDir.append("/") ;
4565  itosConv << msg->grpIndx ;
4566  baseDir += itosConv.str() ;
4567 
4568  if (stat(msg->baseDir, &info) != 0 ) {
4569  CkPrintf( "Node %d cannot access directory %s\n",
4570  CkMyPe(), baseDir.c_str() );
4571  NAMD_die("QM calculation could not be ran. Check your qmBaseDir!");
4572  }
4573  else if (! (stat(baseDir.c_str(), &info) == 0 && S_ISDIR(info.st_mode)) ) {
4574  DebugM(4,"Creating directory " << baseDir.c_str() << std::endl);
4575  int retVal = mkdir(baseDir.c_str(), S_IRWXU);
4576  }
4577 
4578  #ifdef DEBUG_QM
4579  Write_PDB(std::string(baseDir)+"/input.pdb", msg ) ;
4580  #endif
4581 
4582  inputFileName.clear();
4583  inputFileName.append(baseDir.c_str()) ;
4584  inputFileName.append("/qmmm_") ;
4585  inputFileName += itosConv.str() ;
4586  inputFileName.append(".input") ;
4587 
4588  // Opens file for coordinate and parameter input
4589  inputFile = fopen(inputFileName.c_str(),"w");
4590  if ( ! inputFile ) {
4591  iout << iERROR << "Could not open input file for writing: "
4592  << inputFileName << "\n" << endi ;
4593  NAMD_err("Error writing QM input file.");
4594  }
4595 
4596  // Builds the command that will be executed
4597  qmCommand.clear();
4598  qmCommand.append("cd ");
4599  qmCommand.append(baseDir);
4600  qmCommand.append(" ; ");
4601  qmCommand.append(msg->execPath) ;
4602  qmCommand.append(" ") ;
4603  qmCommand.append(inputFileName) ;
4604 
4605  // Builds the file name where orca will place the gradient
4606  // This will be relative to the input file
4607  outputFileName = inputFileName ;
4608  outputFileName.append(".result") ;
4609 
4610  int numPntChrgs = 0;
4611  for (int i=0; i<msg->numAllPntChrgs; i++ ) {
4612  if (pcP[i].type != QMPCTYPE_IGNORE)
4613  numPntChrgs++;
4614  }
4615 
4616  iret = fprintf(inputFile,"%d %d\n",msg->numAllAtoms, numPntChrgs);
4617  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
4618 
4619  DebugM(4, "Writing " << msg->numAllAtoms << " QM atom coords in file " <<
4620  inputFileName.c_str() << " and " << msg->numAllPntChrgs <<
4621  " point charges." << std::endl);
4622 
4623  // write QM and dummy atom coordinates to input file.
4624  QMAtomData *atmP = msg->data ;
4625  for (size_t i=0; i<msg->numAllAtoms; ++i, ++atmP ) {
4626 
4627  double x = atmP->position.x;
4628  double y = atmP->position.y;
4629  double z = atmP->position.z;
4630 
4631  iret = fprintf(inputFile,"%f %f %f %s\n",
4632  x,y,z,atmP->element);
4633  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
4634 
4635  }
4636 
4637  int numWritenPCs = 0;
4638  // Write point charges to file.
4639  pcP = msg->data + msg->numAllAtoms ;
4640  for ( size_t j=0; j < msg->numAllPntChrgs; j++, ++pcP) {
4641 
4642  if (pcP->type == QMPCTYPE_IGNORE)
4643  continue;
4644 
4645  double charge = pcP->charge;
4646 
4647  double x = pcP->position.x;
4648  double y = pcP->position.y;
4649  double z = pcP->position.z;
4650 
4651  iret = fprintf(inputFile,"%f %f %f %f\n",
4652  x,y,z,charge);
4653  if ( iret < 0 ) { NAMD_err("Error writing QM input file."); }
4654 
4655  numWritenPCs++;
4656  }
4657 
4658  DebugM(4,"Closing input file\n");
4659  fclose(inputFile);
4660 
4661  if (msg->prepProcOn) {
4662 
4663  std::string prepProc(msg->prepProc) ;
4664  prepProc.append(" ") ;
4665  prepProc.append(inputFileName) ;
4666  iret = system(prepProc.c_str());
4667  if ( iret == -1 ) { NAMD_err("Error running preparation command for QM calculation."); }
4668  if ( iret ) { NAMD_die("Error running preparation command for QM calculation."); }
4669  }
4670 
4671  // runs QM command
4672  DebugM(4,"Running command ->" << qmCommand.c_str() << "<-" << std::endl);
4673  iret = system(qmCommand.c_str());
4674 
4675  if ( iret == -1 ) { NAMD_err("Error running command for QM forces calculation."); }
4676  if ( iret ) { NAMD_die("Error running command for QM forces calculation."); }
4677 
4678  if (msg->secProcOn) {
4679 
4680  std::string secProc(msg->secProc) ;
4681  secProc.append(" ") ;
4682  secProc.append(inputFileName) ;
4683  itosConv.str("");
4684  itosConv.clear() ;
4685  itosConv << msg->timestep ;
4686  secProc.append(" ") ;
4687  secProc += itosConv.str() ;
4688 
4689  iret = system(secProc.c_str());
4690  if ( iret == -1 ) { NAMD_err("Error running second command for QM calculation."); }
4691  if ( iret ) { NAMD_die("Error running second command for QM calculation."); }
4692  }
4693 
4694  // remove coordinate file
4695 // iret = remove(inputFileName);
4696 // if ( iret ) { NAMD_err(0); }
4697 
4698  // remove coordinate file
4699 // iret = remove(pntChrgFileName);
4700 // if ( iret ) { NAMD_err(0); }
4701 
4702  // opens output file
4703  DebugM(4,"Reading QM data from file " << outputFileName.c_str() << std::endl);
4704  outputFile = fopen(outputFileName.c_str(),"r");
4705  if ( ! outputFile ) {
4706  iout << iERROR << "Could not find QM output file!\n" << endi;
4707  NAMD_err(0);
4708  }
4709 
4710  // Resets the pointers.
4711  atmP = msg->data ;
4712  pcP = msg->data + msg->numAllAtoms ;
4713 
4714  // Allocates the return message.
4715  QMGrpResMsg *resMsg = new (msg->numQMAtoms + msg->numRealPntChrgs, 0) QMGrpResMsg;
4716  resMsg->grpIndx = msg->grpIndx;
4717  resMsg->numForces = msg->numQMAtoms + msg->numRealPntChrgs;
4718  resMsg->energyOrig = 0;
4719  resMsg->energyCorr = 0;
4720  for ( int k=0; k<3; ++k )
4721  for ( int l=0; l<3; ++l )
4722  resMsg->virial[k][l] = 0;
4723  QMForce *resForce = resMsg->force ;
4724 
4725  // We split the data into two pointers so we can "skip" the dummy atoms
4726  // which have no representation in NAMD.
4727  for (int i=0; i<resMsg->numForces; i++) {
4728  resForce[i].force = 0;
4729  resForce[i].charge = 0 ;
4730  if (i < msg->numQMAtoms) {
4731  // We use the replace field to indicate QM atoms,
4732  // which will have charge information.
4733  resForce[i].replace = 1 ;
4734  resForce[i].id = atmP->id;
4735  atmP++;
4736  }
4737  else {
4738  // We use the replace field to indicate QM atoms,
4739  // which will have charge information.
4740  resForce[i].replace = 0 ;
4741  resForce[i].id = pcP->id;
4742  pcP++;
4743  }
4744  }
4745 
4746  // Resets the pointers.
4747  atmP = msg->data ;
4748  pcP = msg->data + msg->numAllAtoms ;
4749 
4750  // Reads the data form the output file created by the QM software.
4751  // Gradients over the QM atoms, and Charges for QM atoms will be read.
4752 
4753  // Number of point charges for which we will receive forces.
4754  int usrPCnum = 0;
4755  const size_t lineLen = 256;
4756  char *line = new char[lineLen];
4757 
4758  fgets(line, lineLen, outputFile);
4759 
4760 // iret = fscanf(outputFile,"%lf %d\n", &resMsg->energyOrig, &usrPCnum);
4761  iret = sscanf(line,"%lf %i\n", &resMsg->energyOrig, &usrPCnum);
4762  if ( iret < 1 ) {
4763  NAMD_die("Error reading energy from QM results file.");
4764  }
4765 
4766  resMsg->energyCorr = resMsg->energyOrig;
4767 
4768  if (iret == 2 && numWritenPCs != usrPCnum) {
4769  iout << iERROR << "Number of point charges does not match what was provided!\n" << endi ;
4770  NAMD_die("Error reading QM results file.");
4771  }
4772 
4773  size_t atmIndx;
4774  double localForce[3];
4775  double localCharge;
4776  for (atmIndx = 0; atmIndx < msg->numAllAtoms; atmIndx++) {
4777 
4778  iret = fscanf(outputFile,"%lf %lf %lf %lf\n",
4779  localForce+0,
4780  localForce+1,
4781  localForce+2,
4782  &localCharge);
4783  if ( iret != 4 ) {
4784  NAMD_die("Error reading forces and charge from QM results file.");
4785  }
4786 
4787  // If we are reading charges and forces on QM atoms, store
4788  // them directly.
4789  if (atmIndx < msg->numQMAtoms ) {
4790 
4791  resForce[atmIndx].force.x = localForce[0];
4792  resForce[atmIndx].force.y = localForce[1];
4793  resForce[atmIndx].force.z = localForce[2];
4794 
4795  atmP[atmIndx].charge = localCharge;
4796  resForce[atmIndx].charge = localCharge;
4797  }
4798 
4799  // If we are reading charges and forces applied on Dummy atoms,
4800  // place them on the appropriate QM and MM atoms to conserve energy.
4801 
4802  // This implementation was based on the description in
4803  // DOI: 10.1002/jcc.20857 : Equations 30 to 32
4804  if ( msg->numQMAtoms <= atmIndx &&
4805  atmIndx < msg->numAllAtoms ) {
4806 
4807  // We keep the calculated charge in the dummy atom
4808  // so that we can calculate electrostatic forces later on.
4809  atmP[atmIndx].charge = localCharge;
4810 
4811  // If we are reading charges from Dummy atoms,
4812  // place them on the appropriate QM atom.
4813  // The dummy atom points to the QM atom to which it is bound.
4814  int qmInd = atmP[atmIndx].bountToIndx ;
4815  resForce[qmInd].charge += localCharge;
4816 
4817  // The dummy atom points to the QM atom to which it is bound.
4818  // The QM atom and the MM atom (in a QM-MM bond) point to each other.
4819  int mmInd = atmP[qmInd].bountToIndx ;
4820 
4821  Vector dir = pcP[mmInd].position - atmP[qmInd].position ;
4822  Real mmqmDist = dir.length() ;
4823 
4824  Real linkDist = Vector(atmP[atmIndx].position -
4825  atmP[qmInd].position).length() ;
4826 
4827  Force mmForce(0), qmForce(0),
4828  linkForce(localForce[0], localForce[1], localForce[2]);
4829 
4830  Vector base = (linkDist/(mmqmDist*mmqmDist*mmqmDist))*dir;
4831  // Unit vectors
4832  Vector xuv(1,0,0), yuv(0,1,0), zuv(0,0,1);
4833  Real xDelta = pcP[mmInd].position.x - atmP[qmInd].position.x;
4834  Real yDelta = pcP[mmInd].position.y - atmP[qmInd].position.y;
4835  Real zDelta = pcP[mmInd].position.z - atmP[qmInd].position.z;
4836 
4837  qmForce += (linkForce*((1 - linkDist/mmqmDist)*xuv +
4838  (xDelta)*base) )*xuv;
4839 
4840  qmForce += (linkForce*((1 - linkDist/mmqmDist)*yuv +
4841  (yDelta)*base) )*yuv;
4842 
4843  qmForce += (linkForce*((1 - linkDist/mmqmDist)*zuv +
4844  (zDelta)*base) )*zuv;
4845 
4846 
4847  mmForce += (linkForce*((linkDist/mmqmDist)*xuv -
4848  (xDelta)*base) )*xuv;
4849 
4850  mmForce += (linkForce*((linkDist/mmqmDist)*yuv -
4851  (yDelta)*base) )*yuv;
4852 
4853  mmForce += (linkForce*((linkDist/mmqmDist)*zuv -
4854  (zDelta)*base) )*zuv;
4855 
4856  resForce[qmInd].force += qmForce;
4857  resForce[msg->numQMAtoms + mmInd].force += mmForce;
4858  }
4859  }
4860 
4861  // The initial point charge index for the force message is the number of QM
4862  // atoms, since the dummy atoms have no representation in NAMD
4863  int pcIndx = msg->numQMAtoms;
4864 
4865  if (usrPCnum > 0) {
4866  // We loop over point charges, reading the total electrostatic force
4867  // applied on them by the QM region.
4868  // We redistribute the forces applied over virtual point
4869  // charges to the MM1 and MM2 atoms (if any virtual PCs exists).
4870  for (size_t i=0; i < msg->numAllPntChrgs; i++, pcIndx++ ) {
4871 
4872  Force totalForce(0);
4873 
4874  // No force was applied to the QM region due to this charge, since it
4875  // was skipped when writing down point charges to the QM software, so it
4876  // does not receive any back from the QM region. It must be an MM1 atom
4877  // from a QM-MM bond.
4878  if (pcP[i].type == QMPCTYPE_IGNORE)
4879  continue;
4880 
4881  iret = fscanf(outputFile,"%lf %lf %lf\n",
4882  &totalForce[0], &totalForce[1], &totalForce[2]);
4883  if ( iret != 3 ) {
4884  NAMD_die("Error reading PC forces from QM results file.");
4885  }
4886 
4887  if (pcP[i].type == QMPCTYPE_CLASSICAL) {
4888  // Checking pcP was not a QM atom in another region
4889  // If so the interaction is accounted in the other region
4890  if (qmAtomGroup[pcP[i].id] == 0) {
4891  // If we already ignored MM1 charges, we take all other
4892  // non-virtual charges and apply forces directly to them.
4893  resForce[pcIndx].force += totalForce;
4894  }
4895  }
4896  else {
4897  // If we are handling virtual PC, we distribute the force over
4898  // MM1 and MM2.
4899 
4900  Force mm1Force(0), mm2Force(0);
4901 
4902  // Virtual PC are bound to MM2.
4903  int mm2Indx = pcP[i].bountToIndx;
4904  // MM2 charges are bound to MM1.
4905  int mm1Indx = pcP[mm2Indx].bountToIndx;
4906 
4907  Real Cq = pcP[i].dist;
4908 
4909  mm1Force = (1-Cq)*totalForce ;
4910  mm2Force = Cq*totalForce ;
4911 
4912  resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
4913  resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
4914  }
4915 
4916 
4917  }
4918  }
4919 
4920  fclose(outputFile);
4921  delete [] line;
4922 
4923  // In case charges are not to be read form the QM software output,
4924  // we load the origianl atom charges.
4925  if (msg->qmAtmChrgMode == QMCHRGNONE) {
4926  int globalNumQMAtoms = Node::Object()->molecule->get_numQMAtoms();
4927  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
4928  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
4929 
4930  atmIndx = 0 ;
4931  for (; atmIndx < msg->numQMAtoms; atmIndx++) {
4932 
4933  // Loops over all QM atoms (in all QM groups) comparing their global indices
4934  for (int qmIter=0; qmIter<globalNumQMAtoms; qmIter++) {
4935 
4936  if (qmAtmIndx[qmIter] == atmP[atmIndx].id) {
4937 
4938  atmP[atmIndx].charge = qmAtmChrg[qmIter];
4939  resForce[atmIndx].charge = qmAtmChrg[qmIter];
4940 
4941  break;
4942  }
4943 
4944  }
4945 
4946  }
4947  for (size_t j=msg->numQMAtoms; j<msg->numAllAtoms; ++j ) {
4948  atmP[j].charge = 0;
4949  }
4950  }
4951 
4952  // remove force file
4953 // DebugM(4, "Removing output file: " << outputFileName << std::endl) ;
4954 // iret = remove(outputFileName);
4955 // if ( iret ) { NAMD_err(0); }
4956 
4957  if (usrPCnum == 0) {
4958  DebugM(4, "Applying forces on " << msg->numRealPntChrgs << " point charges" << std::endl) ;
4959 
4960  atmP = msg->data ;
4961  pcP = msg->data + msg->numAllAtoms ;
4962 
4963  // We only loop over point charges from real atoms, ignoring the ones
4964  // created to handle QM-MM bonds.
4965  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
4966 
4967  // No force was applied to the QM region due to this charge, so it
4968  // does not receive any back from the QM region. It must be an MM1 atom
4969  // from a QM-MM bond.
4970  if (pcP[i].type == QMPCTYPE_IGNORE)
4971  continue;
4972 
4973  Force totalForce(0);
4974 
4975  BigReal pntCharge = pcP[i].charge;
4976 
4977  Position posMM = pcP[i].position ;
4978 
4979  for (size_t j=0; j<msg->numAllAtoms; ++j ) {
4980 
4981  BigReal qmCharge = atmP[j].charge ;
4982 
4983  BigReal force = pntCharge*qmCharge*constants ;
4984 
4985  Position rVec = posMM - atmP[j].position ;
4986 
4987  force /= rVec.length2();
4988 
4989  // We accumulate the total force felt by a point charge
4990  // due to the QM charge distribution. This total force
4991  // will be applied to the point charge if it is a real one,
4992  // or will be distirbuted over MM1 and MM2 point charges, it
4993  // this is a virtual point charge.
4994  totalForce += force*rVec.unit();
4995  }
4996 
4997  if (pcP[i].type == QMPCTYPE_CLASSICAL) {
4998  // Checking pcP was not a QM atom in another region
4999  // If so the interaction is accounted in the other region
5000  if (qmAtomGroup[pcP[i].id] == 0) {
5001  // If we already ignored MM1 charges, we take all other
5002  // non-virtual charges and apply forces directly to them.
5003  resForce[pcIndx].force += totalForce;
5004  }
5005  }
5006  else {
5007  // If we are handling virtual PC, we distribute the force over
5008  // MM1 and MM2.
5009 
5010  Force mm1Force(0), mm2Force(0);
5011 
5012  // Virtual PC are bound to MM2.
5013  int mm2Indx = pcP[i].bountToIndx;
5014  // MM2 charges are bound to MM1.
5015  int mm1Indx = pcP[mm2Indx].bountToIndx;
5016 
5017  Real Cq = pcP[i].dist;
5018 
5019  mm1Force = (1-Cq)*totalForce ;
5020  mm2Force = Cq*totalForce ;
5021 
5022  resForce[msg->numQMAtoms + mm1Indx].force += mm1Force;
5023  resForce[msg->numQMAtoms + mm2Indx].force += mm2Force;
5024  }
5025 
5026  }
5027  }
5028 
5029  // Adjusts forces from PME, canceling contributions from the QM and
5030  // direct Coulomb forces calculated here.
5031  if (msg->PMEOn) {
5032 
5033  DebugM(1,"Correcting forces and energy for PME.\n");
5034 
5035  ewaldcof = msg->PMEEwaldCoefficient;
5036  BigReal TwoBySqrtPi = 1.12837916709551;
5037  pi_ewaldcof = TwoBySqrtPi * ewaldcof;
5038 
5039  for (size_t i=0; i < msg->numQMAtoms; i++) {
5040 
5041  BigReal p_i_charge = atmP[i].charge ;
5042  Position pos_i = atmP[i].position ;
5043 
5044  for (size_t j=i+1; j < msg->numQMAtoms; j++) {
5045 
5046  BigReal p_j_charge = atmP[j].charge ;
5047 
5048  Position pos_j = atmP[j].position ;
5049 
5050  BigReal r = Vector(pos_i - pos_j).length() ;
5051 
5052  BigReal tmp_a = r * ewaldcof;
5053  BigReal tmp_b = erfc(tmp_a);
5054  BigReal corr_energy = tmp_b;
5055  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
5056 
5057 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
5058  BigReal recip_energy = (1-tmp_b)/r;
5059 
5060 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
5061  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
5062 
5063  const BigReal kq_i = p_i_charge * constants;
5064 
5065  // Final force and energy correction for this pair of atoms.
5066  BigReal energy = kq_i * p_j_charge * recip_energy ;
5067 
5068  Force fixForce = -1*kq_i*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
5069 
5070  // The force is *subtracted* from the total force acting on
5071  // both atoms. The sign on fixForce corrects the orientation
5072  // of the subtracted force.
5073 // DebugM(4,"Old forces for QM " << i << ": " << resForce[i].force
5074 // << std::endl);
5075 // DebugM(4,"Old forces for QM " << j << ": " << resForce[j].force
5076 // << std::endl);
5077 // DebugM(4,"Force correction: " << fixForce << std::endl);
5078  resForce[i].force -= fixForce ;
5079  resForce[j].force -= -1*fixForce;
5080 
5081  // The energy is *subtracted* from the total energy calculated here.
5082  resMsg->energyCorr -= energy;
5083  }
5084  }
5085 
5086  pcIndx = msg->numQMAtoms;
5087  // We only loop over point charges from real atoms, ignoring the ones
5088  // created to handle QM-MM bonds.
5089  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
5090 
5091  BigReal p_i_charge = pcPpme[i].charge;
5092  Position pos_i = pcPpme[i].position ;
5093 
5094  const BigReal kq_i = p_i_charge * constants;
5095 
5096  Force fixForce = 0;
5097 
5098  for (size_t j=0; j<msg->numQMAtoms; ++j ) {
5099 
5100  BigReal p_j_charge = atmP[j].charge ;
5101 
5102  Position pos_j = atmP[j].position ;
5103 
5104  BigReal r = Vector(pos_i - pos_j).length() ;
5105 
5106  BigReal tmp_a = r * ewaldcof;
5107  BigReal tmp_b = erfc(tmp_a);
5108  BigReal corr_energy = tmp_b;
5109  BigReal corr_gradient = pi_ewaldcof*exp(-(tmp_a*tmp_a))*r + tmp_b;
5110 
5111 // BigReal recip_energy = (1-tmp_b)/r = erf(tmp_a)/r;
5112  BigReal recip_energy = (1-tmp_b)/r;
5113 
5114 // BigReal recip_gradient = -(1-corr_gradient)/(r*2);
5115  BigReal recip_gradient = -(1-corr_gradient)/(r*r);
5116 
5117  // Final force and energy correction for this pair of atoms.
5118  BigReal energy = kq_i * p_j_charge * recip_energy ;
5119 
5120  fixForce += -1*p_j_charge*(recip_gradient/r)*(pos_i - pos_j) ;
5121 
5122  // The energy is *subtracted* from the total energy calculated here.
5123  resMsg->energyCorr -= energy;
5124 
5125  }
5126 
5127  // The force is *subtracted* from the total force acting on
5128  // the point charge.
5129 // DebugM(4,"Old forces for PC " << pcIndx << ": " << resForce[pcIndx].force
5130 // << std::endl);
5131 // DebugM(4,"Force correction: " << fixForce << std::endl);
5132  resForce[pcIndx].force -= kq_i*fixForce ;
5133  }
5134 
5135  }
5136 
5137  DebugM(1,"Determining virial...\n");
5138 
5139  // Calculates the virial contribution form the forces on QM atoms and
5140  // point charges calculated here.
5141  for (size_t i=0; i < msg->numQMAtoms; i++) {
5142  // virial used by NAMD is -'ve of normal convention, so reverse it!
5143  // virial[i][j] in file should be sum of -1 * f_i * r_j
5144  for ( int k=0; k<3; ++k )
5145  for ( int l=0; l<3; ++l )
5146  resMsg->virial[k][l] = -1*resForce[i].force[k]*atmP[i].position[l];
5147  }
5148 
5149  pcIndx = msg->numQMAtoms; // Index in the real PC force array.
5150  for (size_t i=0; i < msg->numRealPntChrgs; i++, pcIndx++ ) {
5151  // virial used by NAMD is -'ve of normal convention, so reverse it!
5152  // virial[i][j] in file should be sum of -1 * f_i * r_j
5153  for ( int k=0; k<3; ++k )
5154  for ( int l=0; l<3; ++l )
5155  resMsg->virial[k][l] = -1*resForce[pcIndx].force[k]*pcP[i].position[l];
5156  }
5157 
5158  DebugM(1,"End of QM processing. Sending result message.\n");
5159 
5160  // Send message to rank zero with results.
5161  QMProxy[0].recvQMRes(resMsg);
5162 
5163  if (msg->switching && msg->PMEOn)
5164  delete [] pcPpme;
5165 
5166  delete msg;
5167  return ;
5168 }
Node::Object
static Node * Object()
Definition: Node.h:86
QMAtomData::bountToIndx
int bountToIndx
Definition: ComputeQM.C:299
QMGrpCalcMsg::PMEEwaldCoefficient
BigReal PMEEwaldCoefficient
Definition: ComputeQM.C:336
QMGrpResMsg::force
QMForce * force
Definition: ComputeQM.C:181
QMGrpCalcMsg::baseDir
char baseDir[256]
Definition: ComputeQM.C:338
QMGrpCalcMsg::execPath
char execPath[256]
Definition: ComputeQM.C:338
NAMD_err
void NAMD_err(const char *err_msg)
Definition: common.C:106
QMGrpCalcMsg::data
QMAtomData * data
Definition: ComputeQM.C:339
Molecule::get_qmAtmIndx
const int * get_qmAtmIndx()
Definition: Molecule.h:824
QMGrpCalcMsg::constants
BigReal constants
Definition: ComputeQM.C:328
QMAtomData::charge
float charge
Definition: ComputeQM.C:297
Vector
Definition: Vector.h:64
QMGrpCalcMsg::prepProcOn
bool prepProcOn
Definition: ComputeQM.C:330
Molecule::get_numQMAtoms
int get_numQMAtoms()
Definition: Molecule.h:826
Real
float Real
Definition: common.h:109
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
QMGrpCalcMsg::numAllAtoms
int numAllAtoms
Definition: ComputeQM.C:324
QMPCTYPE_CLASSICAL
#define QMPCTYPE_CLASSICAL
Definition: ComputeQM.C:198
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
Vector::z
BigReal z
Definition: Vector.h:66
QMAtomData::dist
Real dist
Definition: ComputeQM.C:302
QMAtomData::position
Position position
Definition: ComputeQM.C:296
iout
#define iout
Definition: InfoStream.h:51
Vector::length
BigReal length(void) const
Definition: Vector.h:169
QMGrpCalcMsg::switching
bool switching
Definition: ComputeQM.C:332
Molecule::get_qmAtomGroup
const Real * get_qmAtomGroup() const
Definition: Molecule.h:820
QMGrpResMsg::numForces
int numForces
Definition: ComputeQM.C:180
QMGrpResMsg::energyOrig
BigReal energyOrig
Definition: ComputeQM.C:177
QMGrpCalcMsg::numQMAtoms
int numQMAtoms
Definition: ComputeQM.C:323
QMGrpCalcMsg::qmAtmChrgMode
int qmAtmChrgMode
Definition: ComputeQM.C:337
QMPCTYPE_IGNORE
#define QMPCTYPE_IGNORE
Definition: ComputeQM.C:197
z
gridSize z
Definition: ComputePmeCUDAKernel.cu:431
QMGrpCalcMsg::numRealPntChrgs
int numRealPntChrgs
Definition: ComputeQM.C:325
Molecule::get_qmAtmChrg
Real * get_qmAtmChrg()
Definition: Molecule.h:823
QMGrpCalcMsg::prepProc
char prepProc[256]
Definition: ComputeQM.C:338
Vector::x
BigReal x
Definition: Vector.h:66
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:85
Vector::length2
BigReal length2(void) const
Definition: Vector.h:173
QMGrpCalcMsg::secProc
char secProc[256]
Definition: ComputeQM.C:338
Vector::y
BigReal y
Definition: Vector.h:66
charge
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
Definition: ComputeNonbondedBase2.h:179
y
gridSize y
Definition: ComputePmeCUDAKernel.cu:430
iERROR
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:83
QMGrpCalcMsg::numAllPntChrgs
int numAllPntChrgs
Definition: ComputeQM.C:326
x
gridSize x
Definition: ComputePmeCUDAKernel.cu:429
QMGrpResMsg::virial
BigReal virial[3][3]
Definition: ComputeQM.C:179
Node::molecule
Molecule * molecule
Definition: Node.h:176
QMGrpCalcMsg::secProcOn
bool secProcOn
Definition: ComputeQM.C:329
Vector::unit
Vector unit(void) const
Definition: Vector.h:182
BigReal
double BigReal
Definition: common.h:114
QMCHRGNONE
#define QMCHRGNONE
Definition: Molecule.h:133
QMGrpResMsg::energyCorr
BigReal energyCorr
Definition: ComputeQM.C:178
for
for(int i=0;i< n1;++i)
Definition: ComputeFullDirectBase.h:34
SortedArray<LSSSubsDat>& ComputeQMMgr::get_subsArray ( )
inline

Definition at line 428 of file ComputeQM.C.

Referenced by lssSubs().

428 { return subsArray; } ;
void ComputeQMMgr::procQMRes ( )

Definition at line 2485 of file ComputeQM.C.

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

2485  {
2486 
2487  // Writes a DCD file with the charges of all QM atoms at a frequency
2488  // defined by the user in qmOutFreq.
2489  if ( simParams->qmOutFreq > 0 &&
2490  timeStep % simParams->qmOutFreq == 0 ) {
2491 
2492  iout << iINFO << "Writing QM charge output at step "
2493  << timeStep << "\n" << endi;
2494 
2495  Real *x = outputData,
2496  *y = outputData + molPtr->get_numQMAtoms(),
2497  *z = outputData + 2*molPtr->get_numQMAtoms();
2498 
2499  for (int qmIt=0; qmIt<numQMAtms; qmIt++){
2500  x[qmIt] = qmCoord[qmIt].id;
2501  y[qmIt] = force[qmCoord[qmIt].id].charge ;
2502  z[qmIt] = 0;
2503  }
2504 
2505  write_dcdstep(dcdOutFile, numQMAtms, x, y, z, 0) ;
2506  }
2507 
2508  // Writes a DCD file with the positions of all QM atoms at a frequency
2509  // defined by the user in qmPosOutFreq.
2510  if ( simParams->qmPosOutFreq > 0 &&
2511  timeStep % simParams->qmPosOutFreq == 0 ) {
2512 
2513  iout << iINFO << "Writing QM position output at step "
2514  << timeStep << "\n" << endi;
2515 
2516  SortedArray<idIndxStr> idIndx;
2517 
2518  for(int i=0; i<numQMAtms;i++) {
2519  idIndx.insert( idIndxStr(qmCoord[i].id, i) );
2520  }
2521  idIndx.sort();
2522 
2523  Real *x = outputData,
2524  *y = outputData + molPtr->get_numQMAtoms(),
2525  *z = outputData + 2*molPtr->get_numQMAtoms();
2526 
2527  for (int qmIt=0; qmIt<numQMAtms; qmIt++){
2528  x[qmIt] = qmCoord[idIndx[qmIt].indx].position.x;
2529  y[qmIt] = qmCoord[idIndx[qmIt].indx].position.y;
2530  z[qmIt] = qmCoord[idIndx[qmIt].indx].position.z;
2531  }
2532 
2533  write_dcdstep(dcdPosOutFile, numQMAtms, x, y, z, 0) ;
2534  }
2535 
2536  // distribute forces
2537  DebugM(4,"Distributing QM forces for all ranks.\n");
2538  for ( int j=0; j < numSources; ++j ) {
2539 
2540  std::set<int> auxset0;
2541  DebugM(1,"Sending forces and charges to source " << j << std::endl);
2542 
2543  QMCoordMsg *qmmsg = 0;
2544  QMPntChrgMsg *pcmsg = 0 ;
2545 
2546  int totForces = 0;
2547  int sourceNode = -1;
2548 
2549  if (pntChrgCoordMsgs == NULL) {
2550 
2551  qmmsg = qmCoordMsgs[j];
2552  qmCoordMsgs[j] = 0;
2553 
2554  totForces = qmmsg->numAtoms ;
2555 
2556  sourceNode = qmmsg->sourceNode;
2557  }
2558  else {
2559  pcmsg = pntChrgCoordMsgs[j];
2560  pntChrgCoordMsgs[j] = 0;
2561 
2562  sourceNode = pcmsg->sourceNode;
2563 
2564  // Since we receive two different messages from nodes, there is no
2565  // guarantee the two sets of messages will come in the same order.
2566  // Therefore, we match the messages by comaring their sourceNodes.
2567  for (int aux=0; aux<numSources; aux++) {
2568 
2569  if (qmCoordMsgs[aux] == 0)
2570  continue;
2571 
2572  qmmsg = qmCoordMsgs[aux];
2573 
2574  if (qmmsg->sourceNode == sourceNode) {
2575  qmCoordMsgs[aux] = 0;
2576  break;
2577  }
2578  }
2579 
2580  DebugM(1,"Building force mesage for rank "
2581  << pcmsg->sourceNode << std::endl);
2582 
2583  totForces = qmmsg->numAtoms + pcmsg->numAtoms;
2584 
2585  // I want to count number of different atom ids
2586  // which is the size of force array.
2587  // Avoids double counting of forces
2588  for ( int i=0; i < qmmsg->numAtoms; ++i ) {
2589  auxset0.insert(qmmsg->coord[i].id);
2590  }
2591  for ( int i=0; i < pcmsg->numAtoms; ++i ) {
2592  auxset0.insert(pcmsg->coord[i].id);
2593  }
2594  totForces = auxset0.size(); // Fixes same patch different QM regions double counting
2595  }
2596 
2597  QMForceMsg *fmsg = new (totForces, subsArray.size(), 0) QMForceMsg;
2598  fmsg->PMEOn = simParams->PMEOn;
2599  fmsg->numForces = totForces;
2600  fmsg->numLssDat = subsArray.size();
2601 
2602  DebugM(1,"Loading QM forces.\n");
2603 
2604  // This iterator is used in BOTH loops!
2605  int forceIter = 0;
2606  auxset0.clear(); // Need to keep track of forces by id that have been copied to fmsg
2607  for ( int i=0; i < qmmsg->numAtoms; ++i ) {
2608  fmsg->force[forceIter] = force[qmmsg->coord[i].id];
2609  forceIter++;
2610  auxset0.insert(qmmsg->coord[i].id); // This should add each qm atom once
2611  }
2612 
2613  delete qmmsg;
2614 
2615  if (pntChrgCoordMsgs != NULL) {
2616  DebugM(1,"Loading PntChrg forces.\n");
2617 
2618  for ( int i=0; i < pcmsg->numAtoms; ++i ) {
2619  // (QM atoms that are PC or repeated PC) are not copied to fmsg
2620  if (auxset0.find(pcmsg->coord[i].id) == auxset0.end()) {
2621  fmsg->force[forceIter] = force[pcmsg->coord[i].id];
2622  forceIter++;
2623  auxset0.insert(pcmsg->coord[i].id);
2624  }
2625  }
2626 
2627  delete pcmsg;
2628  }
2629 
2630  DebugM(1,"A total of " << forceIter << " forces were loaded." << std::endl);
2631 
2632  for ( int i=0; i < subsArray.size(); ++i ) {
2633  fmsg->lssDat[i] = subsArray[i];
2634  }
2635 
2636  #ifdef DEBUG_QM
2637  if (subsArray.size() > 0)
2638  DebugM(3,"A total of " << subsArray.size() << " LSS data structures were loaded." << std::endl);
2639  #endif
2640 
2641  if ( ! j ) {
2642  fmsg->energy = totalEnergy;
2643  for ( int k=0; k<3; ++k )
2644  for ( int l=0; l<3; ++l )
2645  fmsg->virial[k][l] = totVirial[k][l];
2646  } else {
2647  fmsg->energy = 0;
2648  for ( int k=0; k<3; ++k )
2649  for ( int l=0; l<3; ++l )
2650  fmsg->virial[k][l] = 0;
2651  }
2652 
2653  DebugM(4,"Sending forces...\n");
2654 
2655  QMProxy[sourceNode].recvForce(fmsg);
2656 
2657  }
2658 
2659  DebugM(4,"All forces sent from node zero.\n");
2660 }
iINFO
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:81
QMForceMsg::numLssDat
int numLssDat
Definition: ComputeQM.C:190
QMForceMsg::energy
BigReal energy
Definition: ComputeQM.C:187
QMPntChrgMsg::sourceNode
int sourceNode
Definition: ComputeQM.C:169
Molecule::get_numQMAtoms
int get_numQMAtoms()
Definition: Molecule.h:826
QMForceMsg::numForces
int numForces
Definition: ComputeQM.C:189
Real
float Real
Definition: common.h:109
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
Vector::z
BigReal z
Definition: Vector.h:66
iout
#define iout
Definition: InfoStream.h:51
QMCoordMsg::coord
ComputeQMAtom * coord
Definition: ComputeQM.C:110
SortedArray::insert
int insert(const Elem &elem)
SortedArray::sort
void sort(void)
Definition: SortedArray.h:66
ComputeQMAtom::position
Position position
Definition: ComputeQM.C:66
z
gridSize z
Definition: ComputePmeCUDAKernel.cu:431
write_dcdstep
int write_dcdstep(int fd, int N, float *X, float *Y, float *Z, double *cell)
Definition: dcdlib.C:736
QMCoordMsg::sourceNode
int sourceNode
Definition: ComputeQM.C:105
QMForce::charge
float charge
Definition: ComputeQM.h:105
QMForceMsg::force
QMForce * force
Definition: ComputeQM.C:191
Vector::x
BigReal x
Definition: Vector.h:66
QMPntChrgMsg::coord
ComputeQMPntChrg * coord
Definition: ComputeQM.C:171
QMForceMsg::lssDat
LSSSubsDat * lssDat
Definition: ComputeQM.C:192
simParams
#define simParams
Definition: Output.C:127
Vector::y
BigReal y
Definition: Vector.h:66
QMForceMsg::virial
BigReal virial[3][3]
Definition: ComputeQM.C:188
QMCoordMsg::numAtoms
int numAtoms
Definition: ComputeQM.C:106
y
gridSize y
Definition: ComputePmeCUDAKernel.cu:430
ResizeArray::size
int size(void) const
Definition: ResizeArray.h:127
x
gridSize x
Definition: ComputePmeCUDAKernel.cu:429
QMForceMsg::PMEOn
bool PMEOn
Definition: ComputeQM.C:186
void ComputeQMMgr::recvForce ( QMForceMsg *  fmsg)

Definition at line 2662 of file ComputeQM.C.

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

2662  {
2663 
2664  if (CkMyPe()) {
2665  for (int i=0; i<fmsg->numLssDat; i++) {
2666  subsArray.add(fmsg->lssDat[i]) ;
2667  }
2668  }
2669 
2670  QMCompute->saveResults(fmsg);
2671 }
QMForceMsg::numLssDat
int numLssDat
Definition: ComputeQM.C:190
SortedArray::add
int add(const Elem &elem)
Definition: SortedArray.h:55
QMForceMsg::lssDat
LSSSubsDat * lssDat
Definition: ComputeQM.C:192
ComputeQM::saveResults
void saveResults(QMForceMsg *)
Definition: ComputeQM.C:2673
void ComputeQMMgr::recvFullQM ( QMCoordMsg *  qmFullMsg)

Definition at line 1281 of file ComputeQM.C.

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

1281  {
1282 
1283  if (subsArray.size() > 0)
1284  subsArray.clear();
1285 
1286  QMCompute->processFullQM(qmFullMsg);
1287 }
ResizeArray::clear
void clear()
Definition: ResizeArray.h:87
ComputeQM::processFullQM
void processFullQM(QMCoordMsg *)
Definition: ComputeQM.C:1292
ResizeArray::size
int size(void) const
Definition: ResizeArray.h:127
void ComputeQMMgr::recvPartQM ( QMCoordMsg *  msg)

Definition at line 818 of file ComputeQM.C.

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

819 {
820  // In the first (ever) step of the simulation, we allocate arrays that
821  // will be used throughout the simulation, and get some important numbers.
822  if ( ! numSources ) {
823  DebugM(4,"Initializing ComputeQMMgr variables." << std::endl);
824  numSources = (PatchMap::Object())->numNodesWithPatches();
825 
826  DebugM(4,"There are " << numSources << " nodes with patches." << std::endl);
827  qmCoordMsgs = new QMCoordMsg*[numSources];
828  for ( int i=0; i<numSources; ++i ) {
829  qmCoordMsgs[i] = 0;
830  }
831 
832  // Prepares the allocation for the recvPntChrg function.
833 
834  DebugM(4,"Getting data from molecule and simParameters." << std::endl);
835 
836  molPtr = Node::Object()->molecule;
837  simParams = Node::Object()->simParameters;
838 
839  numAtoms = molPtr->numAtoms;
840  force = new QMForce[numAtoms];
841 
842  numQMAtms = molPtr->get_numQMAtoms();
843  qmAtmIter = 0;
844 
845  noPC = simParams->qmNoPC ;
846  meNumMMIndx = molPtr->get_qmMeNumBonds();
847  if (noPC && meNumMMIndx == 0) {
848  pntChrgCoordMsgs = NULL;
849  }
850  else {
851  pntChrgCoordMsgs = new QMPntChrgMsg*[numSources];
852  for ( int i=0; i<numSources; ++i ) {
853  pntChrgCoordMsgs[i] = 0;
854  }
855  }
856 
857  qmPCFreq = molPtr->get_qmPCFreq();
858  resendPCList = false;
859 
860  grpID = molPtr->get_qmGrpID() ;
861  bondValType = simParams->qmBondValType;
862 
863  numQMGrps = molPtr->get_qmNumGrps();
864 
865  grpChrg = molPtr->get_qmGrpChrg() ;
866 
867  grpMult = molPtr->get_qmGrpMult() ;
868 
869  qmLSSOn = simParams->qmLSSOn ;
870  if (qmLSSOn) {
871  qmLSSFreq = molPtr->get_qmLSSFreq() ;
872  qmLSSSize = molPtr->get_qmLSSSize() ;
873  qmLSSIdxs = molPtr->get_qmLSSIdxs() ;
874  qmLSSMass = molPtr->get_qmLSSMass() ;
875  qmLSSResSize = molPtr->get_qmLSSResSize() ;
876  qmLSSRefIDs = molPtr->get_qmLSSRefIDs() ;
877  qmLSSRefSize = molPtr->get_qmLSSRefSize() ;
878 
879  lssPrepare();
880  }
881 
882  numArrivedQMMsg = 0 ;
883  numArrivedPntChrgMsg = 0 ;
884 
885  qmCoord = new ComputeQMAtom[numQMAtms];
886 
887  replaceForces = 0;
888  if (molPtr->get_qmReplaceAll()) {
889  replaceForces = 1;
890  }
891 
892  cSMDon = molPtr->get_qmcSMD() ;
893  if (cSMDon) {
894 
895  // We have to initialize the guide particles during the first step.
896  cSMDInitGuides = 0;
897 
898  cSMDnumInstances = molPtr->get_cSMDnumInst();
899  cSMDindex = molPtr->get_cSMDindex();
900  cSMDindxLen = molPtr->get_cSMDindxLen();
901  cSMDpairs = molPtr->get_cSMDpairs();
902  cSMDKs = molPtr->get_cSMDKs();
903  cSMDVels = molPtr->get_cSMDVels();
904  cSMDcoffs = molPtr->get_cSMDcoffs();
905 
906  cSMDguides = new Position[cSMDnumInstances];
907  cSMDForces = new Force[cSMDnumInstances];
908  }
909 
910 
911  DebugM(4,"Initializing DCD file for charge information." << std::endl);
912 
913  // Initializes output DCD file for charge information.
914  if (simParams->qmOutFreq > 0) {
915  std::string dcdOutFileStr;
916  dcdOutFileStr.clear();
917  dcdOutFileStr.append(simParams->outputFilename) ;
918  dcdOutFileStr.append(".qdcd") ;
919  dcdOutFile = open_dcd_write(dcdOutFileStr.c_str()) ;
920 
921  if (dcdOutFile == DCD_FILEEXISTS) {
922  iout << iERROR << "DCD file " << dcdOutFile << " already exists!!\n" << endi;
923  NAMD_err("Could not write QM charge DCD file.");
924  } else if (dcdOutFile < 0) {
925  iout << iERROR << "Couldn't open DCD file " << dcdOutFile << ".\n" << endi;
926  NAMD_err("Could not write QM charge DCD file.");
927  } else if (! dcdOutFile) {
928  NAMD_bug("ComputeQMMgr::recvPartQM open_dcd_write returned fileid of zero");
929  }
930 
931  timeStep = simParams->firstTimestep;
932 
933  int NSAVC, NFILE, NPRIV, NSTEP;
934  NSAVC = simParams->qmOutFreq;
935  NPRIV = timeStep;
936  NSTEP = NPRIV - NSAVC;
937  NFILE = 0;
938 
939  // Write out the header
940  int ret_code = write_dcdheader(dcdOutFile, dcdOutFileStr.c_str(),
941  numQMAtms, NFILE, NPRIV, NSAVC, NSTEP,
942  simParams->dt/TIMEFACTOR, 0);
943 
944  if (ret_code<0) {
945  NAMD_err("Writing of DCD header failed!!");
946  }
947 
948  // The DCD write function takes 3 independent arrays for X, Y and Z
949  // coordinates, but we allocate one and send in the pieces.
950  outputData = new Real[3*numQMAtms];
951  }
952 
953  DebugM(4,"Initializing DCD file for position information." << std::endl);
954  // Initializes output DCD file for position information.
955  if (simParams->qmPosOutFreq > 0) {
956  std::string dcdPosOutFileStr;
957  dcdPosOutFileStr.clear();
958  dcdPosOutFileStr.append(simParams->outputFilename) ;
959  dcdPosOutFileStr.append(".QMonly.dcd") ;
960  dcdPosOutFile = open_dcd_write(dcdPosOutFileStr.c_str()) ;
961 
962  if (dcdPosOutFile == DCD_FILEEXISTS) {
963  iout << iERROR << "DCD file " << dcdPosOutFile << " already exists!!\n" << endi;
964  NAMD_err("Could not write QM charge DCD file.");
965  } else if (dcdPosOutFile < 0) {
966  iout << iERROR << "Couldn't open DCD file " << dcdPosOutFile << ".\n" << endi;
967  NAMD_err("Could not write QM charge DCD file.");
968  } else if (! dcdPosOutFile) {
969  NAMD_bug("ComputeQMMgr::recvPartQM open_dcd_write returned fileid of zero");
970  }
971 
972  timeStep = simParams->firstTimestep;
973 
974  int NSAVC, NFILE, NPRIV, NSTEP;
975  NSAVC = simParams->qmOutFreq;
976  NPRIV = timeStep;
977  NSTEP = NPRIV - NSAVC;
978  NFILE = 0;
979 
980  // Write out the header
981  int ret_code = write_dcdheader(dcdPosOutFile, dcdPosOutFileStr.c_str(),
982  numQMAtms, NFILE, NPRIV, NSAVC, NSTEP,
983  simParams->dt/TIMEFACTOR, 0);
984 
985  if (ret_code<0) {
986  NAMD_err("Writing of DCD header failed!!");
987  }
988 
989  // The DCD write function takes 3 independent arrays for X, Y and Z
990  // coordinates, but we allocate one and send in the pieces.
991  outputData = new Real[3*numQMAtms];
992  }
993 
994  // Prepares list of PEs which will run the QM software
995  int simsPerNode = simParams->qmSimsPerNode ;
996  int zeroNodeSize = CmiNumPesOnPhysicalNode(0);
997 
998  // Check if the node has enought PEs to run the requested number of simulations.
999  if ( zeroNodeSize < simsPerNode ) {
1000  iout << iERROR << "There are " << zeroNodeSize << " cores pernode, but "
1001  << simsPerNode << " QM simulations per node were requested.\n" << endi ;
1002  NAMD_die("Error preparing QM simulations.");
1003  }
1004 
1005  DebugM(4,"Preparing PE list for QM execution.\n");
1006  qmPEs.clear(); // Making sure its empty.
1007 
1008  int numNodes = CmiNumPhysicalNodes();
1009  int numPlacedQMGrps = 0;
1010  int placedOnNode = 0;
1011  int nodeIt = 0 ;
1012 
1013  // The default is to only run on rank zero.
1014  if ( simsPerNode <= 0 ) {
1015  qmPEs.push_back(0);
1016  numPlacedQMGrps = 1;
1017  }
1018 
1019  while ( (numPlacedQMGrps < numQMGrps) && (simsPerNode > 0) ) {
1020 
1021  // If we searched all nodes, break the loop.
1022  if (nodeIt == numNodes) {
1023  break;
1024  }
1025 
1026  int *pelist;
1027  int nodeSize;
1028  CmiGetPesOnPhysicalNode(CmiPhysicalNodeID(nodeIt), &pelist, &nodeSize);
1029 
1030  DebugM(4,"Checking node " << nodeIt +1 << " out of " << numNodes
1031  << " (" << nodeSize << " PEs: " << pelist[0] << " to "
1032  << pelist[nodeSize-1] << ")." << std::endl );
1033 
1034  for ( int i=0; i<nodeSize; ++i ) {
1035 
1036  // Check if the PE has patches. We only run on PEs with patches!
1037  if ( (PatchMap::Object())->numPatchesOnNode(pelist[i]) == 0 ) {
1038  DebugM(1,"PE " << pelist[i] << " has no patches! We'll skip it."
1039  << std::endl);
1040  continue ;
1041  }
1042 
1043  // Add the target PE on the target node to the list
1044  // of PEs which will carry QM simulations.
1045  qmPEs.push_back(pelist[i]);
1046 
1047  DebugM(1,"Added PE to QM execution list: " << pelist[i] << "\n");
1048 
1049  numPlacedQMGrps++;
1050  placedOnNode++;
1051 
1052  if (placedOnNode == simsPerNode) {
1053  DebugM(1,"Placed enought simulations on this node.\n");
1054  break;
1055  }
1056 
1057 
1058  }
1059 
1060  nodeIt++ ;
1061  placedOnNode = 0;
1062  }
1063 
1064  if ( numPlacedQMGrps < numQMGrps ) {
1065  iout << iWARN << "Could not compute all QM groups in parallel.\n" << endi ;
1066  }
1067 
1068  iout << iINFO << "List of ranks running QM simulations: " << qmPEs[0] ;
1069  for (int i=1; i < qmPEs.size(); i++) {
1070  iout << ", " << qmPEs[i] ;
1071  }
1072  iout << ".\n" << endi;
1073 
1074  }
1075 
1076  DebugM(1,"Receiving from rank " << msg->sourceNode
1077  << " a total of " << msg->numAtoms << " QM atoms." << std::endl);
1078 
1079  // In case we are NOT using point charges but there are QM-MM bonds,
1080  // test each QM message for MM1 atoms.
1081  if (meNumMMIndx > 0) {
1082 
1083  ResizeArray< ComputeQMAtom > tmpAtm;
1084  ComputeQMAtom *data_ptr = msg->coord;
1085 
1086  for (int i=0; i<msg->numAtoms; i++) {
1087  if (data_ptr[i].vdwType < 0) {
1088  tmpAtm.add(data_ptr[i]) ;
1089  }
1090  }
1091 
1092  QMPntChrgMsg *pntChrgMsg = new (tmpAtm.size(), 0) QMPntChrgMsg;
1093  pntChrgMsg->sourceNode = msg->sourceNode ;
1094  pntChrgMsg->numAtoms = tmpAtm.size() ;
1095  ComputeQMPntChrg* newPCData = pntChrgMsg->coord ;
1096 
1097  QMCoordMsg *newMsg = msg;
1098 
1099  if (tmpAtm.size() > 0) {
1100 
1101  newMsg = new (msg->numAtoms - tmpAtm.size(),0, 0) QMCoordMsg;
1102  newMsg->sourceNode = msg->sourceNode ;
1103  newMsg->numAtoms = msg->numAtoms - tmpAtm.size() ;
1104  newMsg->timestep = msg->timestep ;
1105  ComputeQMAtom *newMsgData = newMsg->coord;
1106 
1107  for (int i=0; i<msg->numAtoms; i++) {
1108  if (data_ptr[i].vdwType < 0) {
1109  newPCData->position = data_ptr[i].position ;
1110  newPCData->charge = data_ptr[i].charge ;
1111  newPCData->id = data_ptr[i].id ;
1112  newPCData->qmGrpID = data_ptr[i].qmGrpID ;
1113  newPCData->homeIndx = data_ptr[i].homeIndx ;
1114  newPCData->dist = 0 ;
1115  newPCData->mass = 0 ;
1116  newPCData->vdwType = 0 ;
1117  newPCData++;
1118  }
1119  else {
1120  *newMsgData = data_ptr[i] ;
1121  newMsgData++;
1122  }
1123  }
1124 
1125  delete msg;
1126 
1127  }
1128 
1129  qmCoordMsgs[numArrivedQMMsg] = newMsg;
1130  ++numArrivedQMMsg;
1131 
1132  pntChrgCoordMsgs[numArrivedPntChrgMsg] = pntChrgMsg;
1133  ++numArrivedPntChrgMsg;
1134  }
1135  else {
1136  qmCoordMsgs[numArrivedQMMsg] = msg;
1137  ++numArrivedQMMsg;
1138  }
1139 
1140  if ( numArrivedQMMsg < numSources )
1141  return;
1142 
1143  // Now that all messages arrived, get all QM positions.
1144  for (int msgIt=0; msgIt < numArrivedQMMsg; msgIt++){
1145 
1146  DebugM(1, "Getting positions for " << qmCoordMsgs[msgIt]->numAtoms
1147  << " QM atoms in this message." << std::endl);
1148 
1149  for ( int i=0; i < qmCoordMsgs[msgIt]->numAtoms; ++i ) {
1150  qmCoord[qmAtmIter] = qmCoordMsgs[msgIt]->coord[i];
1151  qmAtmIter++;
1152  }
1153  }
1154 
1155  if (qmAtmIter != numQMAtms) {
1156  iout << iERROR << "The number of QM atoms received (" << qmAtmIter
1157  << ") is different than expected: " << numQMAtms << "\n" << endi;
1158  NAMD_err("Problems broadcasting QM atoms.");
1159  }
1160 
1161  // Resets the counter for the next step.
1162  numArrivedQMMsg = 0;
1163  qmAtmIter = 0;
1164 
1165  timeStep = qmCoordMsgs[0]->timestep;
1166 
1167  // Makes sure there is no trash or old info in the force array.
1168  // This is where we accumulate forces from the QM software and our own
1169  // Coulomb forces. It will have info on QM atoms and point charges only.
1170  for (int i=0; i<numAtoms; ++i ) {
1171  force[i].force = 0; // this assigns 0 (ZERO) to all componenets of the vector.
1172  force[i].replace = 0; // By default, no atom has all forces replaced.
1173  force[i].homeIndx = -1; // To prevent errors from sliping.
1174  force[i].charge = 0;
1175  force[i].id = i; // Initializes the variable since not all forces are sent to all ranks.
1176  }
1177 
1178 
1179  for (int i=0; i<numQMAtms; i++) {
1180  // Each force receives the home index of its atom with respect to the
1181  // local set of atoms in each node.
1182  if (force[qmCoord[i].id].homeIndx != -1
1183  && force[qmCoord[i].id].homeIndx != qmCoord[i].homeIndx
1184  ) {
1185  iout << iERROR << "Overloading QM atom "
1186  << qmCoord[i].id << "; home index: "
1187  << force[qmCoord[i].id].homeIndx << " with " << qmCoord[i].homeIndx
1188  << "\n" << endi ;
1189  NAMD_die("Error preparing QM simulations.");
1190  }
1191 
1192  force[qmCoord[i].id].homeIndx = qmCoord[i].homeIndx;
1193  // Each force on QM atoms has an indicator so NAMD knows if all
1194  // NAMD forces should be erased and completely overwritten by the
1195  // external QM forces.
1196  force[qmCoord[i].id].replace = replaceForces;
1197  }
1198 
1199  if (noPC) {
1200  // this pointer should be freed in rank zero, after receiving it.
1201  QMPntChrgMsg *pntChrgMsg = new (0, 0) QMPntChrgMsg;
1202  pntChrgMsg->sourceNode = CkMyPe();
1203  pntChrgMsg->numAtoms = 0;
1204 
1205  QMProxy[0].recvPntChrg(pntChrgMsg);
1206  }
1207  else {
1208  // The default mode is to update the poitn charge selection at every step.
1209  int pcSelMode = PCMODEUPDATESEL;
1210 
1211  int msgPCListSize = 0;
1212  // We check wether point charges are to be selected with a stride.
1213  if ( qmPCFreq > 0 ) {
1214  if (timeStep % qmPCFreq == 0) {
1215  // Marks that the PC list determined in this step will
1216  // be broadcast on the *next* step, and kept for the following
1217  // qmPCFreq-1 steps.
1218  resendPCList = true;
1219 
1220  // Clears the list since we don't know how many charges
1221  // will be selected this time.
1222  delete [] pcIDList;
1223  }
1224  else {
1225  // If the PC selection is not to be updated in this step,
1226  // inform the nodes that the previous list (or the updated
1227  // list being sent in this message) is to be used and only
1228  // updated positions will be returned.
1229  pcSelMode = PCMODEUPDATEPOS;
1230 
1231  // If this is the first step after a PC re-selection, all
1232  // ranks receive the total list, since atoms may move between
1233  // PEs in between PC re-assignments (if they are far enought apart
1234  // or if you are unlucky)
1235  if (resendPCList) {
1236  msgPCListSize = pcIDListSize;
1237  resendPCList = false;
1238  }
1239  }
1240  }
1241 
1242  // In case we are using custom selection of point charges, indicate the mode.
1243  if (simParams->qmCustomPCSel)
1244  pcSelMode = PCMODECUSTOMSEL;
1245 
1246  DebugM(1,"Broadcasting current positions of QM atoms.\n");
1247  for ( int j=0; j < numSources; ++j ) {
1248  // This pointer will be freed in the receiving rank.
1249  QMCoordMsg *qmFullMsg = new (numQMAtms, msgPCListSize, 0) QMCoordMsg;
1250  qmFullMsg->sourceNode = CkMyPe();
1251  qmFullMsg->numAtoms = numQMAtms;
1252  qmFullMsg->pcSelMode = pcSelMode;
1253  qmFullMsg->numPCIndxs = msgPCListSize;
1254  ComputeQMAtom *data_ptr = qmFullMsg->coord;
1255  int *msgPCListP = qmFullMsg->pcIndxs;
1256 
1257  for (int i=0; i<numQMAtms; i++) {
1258  data_ptr->position = qmCoord[i].position;
1259  data_ptr->charge = qmCoord[i].charge;
1260  data_ptr->id = qmCoord[i].id;
1261  data_ptr->qmGrpID = qmCoord[i].qmGrpID;
1262  data_ptr->homeIndx = -1; // So there is no mistake that there is no info here.
1263  data_ptr++;
1264  }
1265 
1266  for (int i=0; i<msgPCListSize; i++) {
1267  msgPCListP[i] = pcIDList[i] ;
1268  }
1269 
1270  #ifdef DEBUG_QM
1271  if (j == 0)
1272  Write_PDB("qmMsg.pdb", qmFullMsg) ;
1273  #endif
1274 
1275  // The messages are deleted later, we will need them.
1276  QMProxy[qmCoordMsgs[j]->sourceNode].recvFullQM(qmFullMsg);
1277  }
1278  }
1279 }
Node::Object
static Node * Object()
Definition: Node.h:86
QMForce::replace
int replace
Definition: ComputeQM.h:102
Molecule::get_qmLSSRefSize
int * get_qmLSSRefSize()
Definition: Molecule.h:851
iINFO
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:81
QMCoordMsg::timestep
int timestep
Definition: ComputeQM.C:107
NAMD_err
void NAMD_err(const char *err_msg)
Definition: common.C:106
Molecule::get_qmGrpMult
Real * get_qmGrpMult()
Definition: Molecule.h:832
Molecule::get_qmNumGrps
int get_qmNumGrps() const
Definition: Molecule.h:828
PCMODEUPDATEPOS
#define PCMODEUPDATEPOS
Definition: ComputeQM.C:100
PatchMap::Object
static PatchMap * Object()
Definition: PatchMap.h:27
Molecule::get_cSMDnumInst
int get_cSMDnumInst()
Definition: Molecule.h:869
QMPntChrgMsg::sourceNode
int sourceNode
Definition: ComputeQM.C:169
Vector
Definition: Vector.h:64
Node::simParameters
SimParameters * simParameters
Definition: Node.h:178
Molecule::get_numQMAtoms
int get_numQMAtoms()
Definition: Molecule.h:826
Molecule::get_qmLSSMass
Mass * get_qmLSSMass()
Definition: Molecule.h:847
ComputeQMAtom::charge
float charge
Definition: ComputeQM.C:67
QMCoordMsg::pcIndxs
int * pcIndxs
Definition: ComputeQM.C:111
Real
float Real
Definition: common.h:109
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
iWARN
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:82
iout
#define iout
Definition: InfoStream.h:51
Molecule::get_cSMDindxLen
int const * get_cSMDindxLen()
Definition: Molecule.h:870
Molecule::get_qmcSMD
Bool get_qmcSMD()
Definition: Molecule.h:868
Molecule::get_qmLSSRefIDs
int * get_qmLSSRefIDs()
Definition: Molecule.h:850
Molecule::get_qmLSSIdxs
int * get_qmLSSIdxs()
Definition: Molecule.h:846
ComputeQMPntChrg::position
Position position
Definition: ComputeQM.C:126
open_dcd_write
int open_dcd_write(const char *dcdname)
Definition: dcdlib.C:662
Molecule::get_qmGrpChrg
Real * get_qmGrpChrg()
Definition: Molecule.h:831
QMCoordMsg::coord
ComputeQMAtom * coord
Definition: ComputeQM.C:110
Molecule::get_cSMDpairs
int const *const * get_cSMDpairs()
Definition: Molecule.h:872
DCD_FILEEXISTS
#define DCD_FILEEXISTS
Definition: dcdlib.h:52
Molecule::get_qmMeNumBonds
int get_qmMeNumBonds()
Definition: Molecule.h:858
Molecule::get_qmGrpID
Real * get_qmGrpID()
Definition: Molecule.h:830
QMForce::id
int id
Definition: ComputeQM.h:106
PCMODEUPDATESEL
#define PCMODEUPDATESEL
Definition: ComputeQM.C:99
ComputeQMAtom::position
Position position
Definition: ComputeQM.C:66
NAMD_bug
void NAMD_bug(const char *err_msg)
Definition: common.C:129
Molecule::get_qmPCFreq
int get_qmPCFreq()
Definition: Molecule.h:862
Molecule::get_qmLSSFreq
int get_qmLSSFreq()
Definition: Molecule.h:848
QMCoordMsg::sourceNode
int sourceNode
Definition: ComputeQM.C:105
QMForce::charge
float charge
Definition: ComputeQM.h:105
ComputeQMAtom::qmGrpID
Real qmGrpID
Definition: ComputeQM.C:69
Molecule::get_cSMDcoffs
const Real * get_cSMDcoffs()
Definition: Molecule.h:875
Molecule::numAtoms
int numAtoms
Definition: Molecule.h:557
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:85
QMPntChrgMsg::coord
ComputeQMPntChrg * coord
Definition: ComputeQM.C:171
QMForce::force
Force force
Definition: ComputeQM.h:103
write_dcdheader
int write_dcdheader(int fd, const char *filename, int N, int NFILE, int NPRIV, int NSAVC, int NSTEP, double DELTA, int with_unitcell)
Definition: dcdlib.C:915
ResizeArray::add
int add(const Elem &elem)
Definition: ResizeArray.h:97
QMCoordMsg::pcSelMode
int pcSelMode
Definition: ComputeQM.C:109
Molecule::get_cSMDKs
const Real * get_cSMDKs()
Definition: Molecule.h:873
simParams
#define simParams
Definition: Output.C:127
TIMEFACTOR
#define TIMEFACTOR
Definition: common.h:48
QMCoordMsg::numAtoms
int numAtoms
Definition: ComputeQM.C:106
QMForce::homeIndx
int homeIndx
Definition: ComputeQM.h:104
Molecule::get_cSMDVels
const Real * get_cSMDVels()
Definition: Molecule.h:874
Molecule::get_qmLSSSize
int * get_qmLSSSize()
Definition: Molecule.h:845
ResizeArray::size
int size(void) const
Definition: ResizeArray.h:127
iERROR
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:83
PCMODECUSTOMSEL
#define PCMODECUSTOMSEL
Definition: ComputeQM.C:101
Molecule::get_cSMDindex
int const *const * get_cSMDindex()
Definition: Molecule.h:871
Node::molecule
Molecule * molecule
Definition: Node.h:176
Molecule::get_qmReplaceAll
Bool get_qmReplaceAll()
Definition: Molecule.h:855
QMCoordMsg::numPCIndxs
int numPCIndxs
Definition: ComputeQM.C:108
Molecule::get_qmLSSResSize
int get_qmLSSResSize()
Definition: Molecule.h:849
void ComputeQMMgr::recvPntChrg ( QMPntChrgMsg *  msg)

Definition at line 1910 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, ComputeQMAtom::charge, QMAtomData::charge, QMGrpCalcMsg::charge, ResizeArray< T >::clear(), QMGrpCalcMsg::configLines, Node::configList, QMGrpCalcMsg::constants, COULOMB, custom_ComputeQMAtom_Less(), SimParameters::cutoff, QMGrpCalcMsg::cutoff, StringList::data, QMGrpCalcMsg::data, DebugM, SimParameters::dielectric, QMAtomData::dist, QMAtomData::element, endi(), QMGrpCalcMsg::execPath, SortedArray< Type >::find(), ConfigList::find(), SimParameters::firstTimestep, Molecule::get_qmDummyBondVal(), Molecule::get_qmDummyElement(), Molecule::get_qmElements(), Molecule::get_qmGrpBonds(), Molecule::get_qmGrpNumBonds(), Molecule::get_qmGrpSizes(), Molecule::get_qmMMBond(), Molecule::get_qmMMBondedIndx(), Molecule::get_qmMMChargeTarget(), Molecule::get_qmMMNumTargs(), QMGrpCalcMsg::grpIndx, QMForce::homeIndx, QMAtomData::id, iERROR(), iout, SortedArray< Type >::load(), QMGrpCalcMsg::multiplicity, NAMD_die(), StringList::next, SimParameters::nonbondedScaling, QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMPntChrgMsg::numAtoms, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, Node::Object(), QMGrpCalcMsg::pcScheme, QMGrpCalcMsg::peIter, QMGrpCalcMsg::PMEEwaldCoefficient, SimParameters::PMEEwaldCoefficient, QMGrpCalcMsg::PMEOn, SimParameters::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< T >::size(), SortedArray< Type >::sort(), sort, QMGrpCalcMsg::swdist, QMGrpCalcMsg::switching, SimParameters::switchingDist, QMGrpCalcMsg::switchType, QMGrpCalcMsg::timestep, QMAtomData::type, and Vector::unit().

1910  {
1911 
1912  // All the preparation that used to be in this function was moved to
1913  // recvPartQM, which is called first in rank zero.
1914 
1915  if (noPC) {
1916  // Even if we have QM-MM bonds, the point charge messages
1917  // are handled in recvPartQM
1918  delete msg;
1919  }
1920  else {
1921  pntChrgCoordMsgs[numArrivedPntChrgMsg] = msg;
1922  ++numArrivedPntChrgMsg;
1923 
1924  if ( numArrivedPntChrgMsg < numSources )
1925  return;
1926  }
1927 
1928  // Resets counters for next step.
1929  numRecQMRes = 0;
1930 
1931  totalEnergy = 0;
1932 
1933  for ( int k=0; k<3; ++k )
1934  for ( int l=0; l<3; ++l )
1935  totVirial[k][l] = 0;
1936 
1937  // ALL DATA ARRIVED --- CALCULATE FORCES
1938 
1939  const char *const *const elementArray = molPtr->get_qmElements() ;
1940  const char *const *const dummyElmArray = molPtr->get_qmDummyElement();
1941  const int *const qmGrpSize = molPtr->get_qmGrpSizes();
1942 
1943  const BigReal *const dummyBondVal = molPtr->get_qmDummyBondVal();
1944  const int *const grpNumBonds = molPtr->get_qmGrpNumBonds() ;
1945  const int *const *const qmMMBond = molPtr->get_qmMMBond() ;
1946  const int *const *const qmGrpBonds = molPtr->get_qmGrpBonds() ;
1947  const int *const *const qmMMBondedIndx = molPtr->get_qmMMBondedIndx() ;
1948  const int *const *const chargeTarget = molPtr->get_qmMMChargeTarget() ;
1949  const int *const numTargs = molPtr->get_qmMMNumTargs() ;
1950 
1951  BigReal constants = COULOMB * simParams->nonbondedScaling / (simParams->dielectric) ;
1952  // COULOMB is in kcal*Angs/(mol*e^2)
1953 // BigReal constants = COULOMB ;
1954 
1955  if ( qmPCFreq > 0 ) {
1956  DebugM(4,"Using point charge stride of " << qmPCFreq << "\n")
1957  // In case we are skiping steps between PC selection, store a main list
1958  // in rank zero for future steps. Then we only update positions untill
1959  // we reach a new "qmPCFreq" step, when a new PC selection is made.
1960 
1961  if (timeStep % qmPCFreq == 0) {
1962  DebugM(4,"Loading a new selection of PCs.\n")
1963 
1964  // We only re-set this arrya in a step where a new PC selection is made.
1965  pntChrgMsgVec.clear();
1966  for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
1967  // Accumulates the message point charges into a local vector.
1968  for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
1969  pntChrgMsgVec.push_back(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
1970  }
1971  }
1972 
1973  // This fast array is created to send the entire PC IDs list to all
1974  // PEs in the next step.
1975  pcIDListSize = pntChrgMsgVec.size();
1976  pcIDList = new int[pcIDListSize] ;
1977  for (size_t i=0; i<pcIDListSize; i++) {
1978  pcIDList[i] = pntChrgMsgVec[i].id;
1979 
1980  // Loads home indexes of MM atoms received as point charges.
1981  // Each force receives the home index of its atom with respect to the
1982  // local set of atoms in each node.
1983  force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
1984  }
1985  }
1986  else {
1987  DebugM(4,"Updating position/homeIdex of old PC selection.\n")
1988 
1989  // We build a sorted array so that we can quickly access the unordered
1990  // data we just received, and update positions of the same selection
1991  // of point charges.
1992  pcDataSort.clear();
1993  for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
1994  // Accumulates the message point charges into a local sorted array.
1995  for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
1996  pcDataSort.load(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
1997  }
1998  }
1999  pcDataSort.sort();
2000 
2001  iout << "Loaded new positions in sorted array: " << pcDataSort.size() << "\n" << endi;
2002 
2003  // If we are using a set of point charges that was selected in a
2004  // previous step, we update the PC POSITION ONLY.
2005  for (size_t i=0; i<pntChrgMsgVec.size(); i++) {
2006 
2007  auto pntr = pcDataSort.find(pntChrgMsgVec[i]);
2008 
2009  pntChrgMsgVec[i].position = pntr->position ;
2010  pntChrgMsgVec[i].homeIndx = pntr->homeIndx ;
2011 
2012  // Loads home indexes of MM atoms received as point charges.
2013  // Each force receives the home index of its atom with respect to the
2014  // local set of atoms in each node.
2015  force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
2016  }
2017  }
2018  }
2019  else {
2020  DebugM(4,"Updating PCs at every step.\n")
2021 
2022  pntChrgMsgVec.clear();
2023  for (int pcMsgIt=0; pcMsgIt < numArrivedPntChrgMsg; pcMsgIt++) {
2024  // Accumulates the message point charges into a local vector.
2025  for ( int i=0; i < pntChrgCoordMsgs[pcMsgIt]->numAtoms; ++i ) {
2026  pntChrgMsgVec.push_back(pntChrgCoordMsgs[pcMsgIt]->coord[i]);
2027  }
2028  }
2029 
2030  // Loads home indexes of MM atoms received as point charges.
2031  for (size_t i=0; i<pntChrgMsgVec.size(); i++) {
2032  // Each force receives the home index of its atom with respect to the
2033  // local set of atoms in each node.
2034 
2035  #ifdef DEBUG_QM
2036  if (force[pntChrgMsgVec[i].id].homeIndx != -1
2037  and force[pntChrgMsgVec[i].id].homeIndx != pntChrgMsgVec[i].homeIndx
2038  ) {
2039  iout << iERROR << "Overloading point charge "
2040  << pntChrgMsgVec[i].id << "; home index: "
2041  << force[pntChrgMsgVec[i].id].homeIndx << " with " << pntChrgMsgVec[i].homeIndx
2042  << "\n" << endi ;
2043  NAMD_die("Error preparing QM simulations.");
2044  }
2045  #endif
2046 
2047  force[pntChrgMsgVec[i].id].homeIndx = pntChrgMsgVec[i].homeIndx;
2048  }
2049  }
2050 
2051  // Reset counter for next timestep
2052  numArrivedPntChrgMsg = 0;
2053 
2054  DebugM(4,"A total of " << pntChrgMsgVec.size()
2055  << " point charges arrived." << std::endl);
2056 
2057  DebugM(4,"Starting QM groups processing.\n");
2058 
2059  QMAtmVec grpQMAtmVec;
2060  QMPCVec grpPntChrgVec;
2061 
2062  // Final set of charges, created or not, that is sent to the QM software.
2063  // This set will depend on how QM-MM bonds are processed and presented to the
2064  // QM region.
2065  QMPCVec grpAppldChrgVec;
2066 
2067  // Vector of dummy atoms created to treat QM-MM bonds.
2068  std::vector<dummyData> dummyAtoms ;
2069 
2070  // Initializes the loop for receiving the QM results.
2071  thisProxy[0].recvQMResLoop() ;
2072 
2073  // Iterator for target PE where QM simulations will run.
2074  int peIter = 0;
2075 
2076  for (int grpIter = 0; grpIter < numQMGrps; grpIter++) {
2077 
2078  grpQMAtmVec.clear();
2079  grpPntChrgVec.clear();
2080  grpAppldChrgVec.clear();
2081  dummyAtoms.clear();
2082 
2083  DebugM(4,"Calculating QM group " << grpIter +1
2084  << " (ID: " << grpID[grpIter] << ")." << std::endl);
2085 
2086  DebugM(4,"Compiling Config Lines into one string for message...\n");
2087 
2088  // This will hold a big sting with all configuration lines the user supplied.
2089  int lineIter = 0 ;
2090  std::string configLines ;
2091  StringList *current = Node::Object()->configList->find("QMConfigLine");
2092  for ( ; current; current = current->next ) {
2093  std::string confLineStr(current->data);
2094 
2095  // In case we need to add charges to MOPAC command line.
2096  if (simParams->qmFormat == QMFormatMOPAC && simParams->qmMOPACAddConfigChrg && lineIter == 0) {
2097  std::ostringstream itosConv ;
2098  itosConv << grpChrg[grpIter] ;
2099  confLineStr.append( " CHARGE=" );
2100  confLineStr.append( itosConv.str() );
2101 
2102  }
2103 
2104  configLines.append(confLineStr);
2105  configLines.append("\n");
2106 
2107  lineIter++;
2108  }
2109 
2110  DebugM(4,"Determining point charges...\n");
2111 
2112  Real qmTotalCharge = 0;
2113  // Loads the QM atoms for this group.
2114  for (int qmIt=0; qmIt<numQMAtms; qmIt++){
2115  if (qmCoord[qmIt].qmGrpID == grpID[grpIter]) {
2116  grpQMAtmVec.push_back(qmCoord[qmIt]);
2117  qmTotalCharge += qmCoord[qmIt].charge;
2118  }
2119  }
2120  if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
2121  qmTotalCharge = roundf(qmTotalCharge) ;
2122  }
2123 
2124  // Sorts the vector so that the QM message is always built with the same order of atoms.
2125  std::sort(grpQMAtmVec.begin(), grpQMAtmVec.end(), custom_ComputeQMAtom_Less);
2126 
2127  Real pcTotalCharge = 0;
2128  // Loads the point charges to a local vector for this QM group.
2129  for (auto pntChrgIt = pntChrgMsgVec.begin();
2130  pntChrgIt != pntChrgMsgVec.end(); pntChrgIt++) {
2131  if ((*pntChrgIt).qmGrpID == grpID[grpIter] ) {
2132  grpPntChrgVec.push_back( (*pntChrgIt) );
2133  pcTotalCharge += (*pntChrgIt).charge;
2134  }
2135  }
2136  if ((fabsf(roundf(pcTotalCharge) - pcTotalCharge) <= 0.001f) ) {
2137  pcTotalCharge = roundf(pcTotalCharge) ;
2138  }
2139 
2140  #ifdef DEBUG_QM
2141  if (grpQMAtmVec.size() != qmGrpSize[grpIter]) {
2142  iout << iERROR << "The number of QM atoms received for group " << grpID[grpIter]
2143  << " does not match the expected: " << grpQMAtmVec.size()
2144  << " vs " << qmGrpSize[grpIter] << "\n" << endi ;
2145 
2146  NAMD_die("Error processing QM group.");
2147  }
2148  #endif
2149 
2150  DebugM(1,"Found " << grpPntChrgVec.size() << " point charges for QM group "
2151  << grpIter << " (ID: " << grpID[grpIter] << "; Num QM atoms: "
2152  << grpQMAtmVec.size() << "; Num QM-MM bonds: "
2153  << grpNumBonds[grpIter] << ")" << std::endl);
2154 
2155  DebugM(1,"Total QM charge: " << qmTotalCharge
2156  << "; Total point-charge charge: " << pcTotalCharge << std::endl);
2157 
2158  // If we have a frequency for LSS update, check if we shoudl do it in
2159  // the current time step.
2160  if ( qmLSSFreq > 0 && ((timeStep + 1) % qmLSSFreq == 0 )) {
2161  lssUpdate(grpIter, grpQMAtmVec, grpPntChrgVec);
2162  }
2163 
2164  // This function checks data and treats the charge (and existence) of
2165  // the MM atoms in and around QM-MM bonds. It is only executed in
2166  // electrostatic embeding QM/MM simulations.
2167  if (! noPC ) {
2168  procBonds(grpNumBonds[grpIter], qmGrpBonds[grpIter],
2169  qmMMBondedIndx[grpIter],
2170  chargeTarget, numTargs,
2171  grpPntChrgVec, grpAppldChrgVec) ;
2172  }
2173  else {
2174  grpAppldChrgVec = grpPntChrgVec;
2175  }
2176 
2177  // For future use, we get the pairs of indexes of QM atoms and MM atoms which are
2178  // bound in QM-MM bonds.
2179  std::vector< std::pair<int,int> > qmPCLocalIndxPairs ;
2180 
2181  // Create and position dummy atoms.
2182  Position mmPos(0), qmPos(0);
2183  for (int dummyIt = 0; dummyIt < grpNumBonds[grpIter]; dummyIt++) {
2184 
2185  int qmMMBondIndx = qmGrpBonds[grpIter][dummyIt] ;
2186 
2187  BigReal bondVal = dummyBondVal[qmMMBondIndx] ;
2188 
2189  int mmAtmIndx = qmMMBond[qmMMBondIndx][0] ;
2190  int qmAtmIndx = qmMMBond[qmMMBondIndx][1] ;
2191 
2192  // Sicne we don't know in which patch/node the QM atom is, or the
2193  // order in which they will arrive in rank zero, we have
2194  // no direct index to it.
2195  #ifdef DEBUG_QM
2196  bool missingQM = true, missingMM = true;
2197  #endif
2198  size_t qmIt ;
2199  for (qmIt=0; qmIt<grpQMAtmVec.size(); qmIt++){
2200  if (grpQMAtmVec[qmIt].id == qmAtmIndx) {
2201  qmPos = grpQMAtmVec[qmIt].position;
2202 
2203  #ifdef DEBUG_QM
2204  missingQM = false;
2205  #endif
2206 
2207  break;
2208  }
2209  }
2210  // The same is true about the MM atom to which the QM atom is bound,
2211  // we must look
2212  size_t pcIt;
2213  for (pcIt=0; pcIt < grpPntChrgVec.size(); pcIt++) {
2214  if (grpPntChrgVec[pcIt].id == mmAtmIndx) {
2215  mmPos = grpPntChrgVec[pcIt].position ;
2216 
2217  #ifdef DEBUG_QM
2218  missingMM = false;
2219  #endif
2220 
2221  break;
2222  }
2223  }
2224 
2225  qmPCLocalIndxPairs.push_back(std::pair<int,int>(qmIt, pcIt) );
2226 
2227  #ifdef DEBUG_QM
2228  // Checks if the MM atom was included as a
2229  // point charge due proximity to a QM atom, and if the QM atom arrived.
2230  if ( missingMM or missingQM ) {
2231  // If it wasn't, there is an error somwhere.
2232 
2233  if (missingMM) {
2234  iout << iERROR << "The MM atom " << mmAtmIndx
2235  << " is bound to a QM atom, but it was not selected as a poitn charge."
2236  << " Check your cutoff radius!\n" << endi ;
2237 
2238  NAMD_die("Error in QM-MM bond processing.");
2239  }
2240  if (missingQM) {
2241  iout << iERROR << "The QM atom " << qmAtmIndx
2242  << " is bound to an MM atom, but it was not sent to rank zero for processing."
2243  << " Check your configuration!\n" << endi ;
2244 
2245  NAMD_die("Error in QM-MM bond processing.");
2246  }
2247  }
2248  #endif
2249 
2250  Vector bondVec = mmPos - qmPos ;
2251 
2252  if (bondValType == QMLENTYPE) {
2253  // If a length is defined by the user, or a default len
2254  // is used, we calculate the unit vector for the displacement
2255  // and multiply by the desired length in order
2256  // to get the final dummy atom position relative to the
2257  // QM atom.
2258  bondVec = bondVec.unit() ;
2259  bondVec *= bondVal ;
2260  }
2261  else if (bondValType == QMRATIOTYPE) {
2262  // If distance a ratio was defined by the user, then
2263  // the displacement vector is multiplied by that ratio
2264  // to get the final dummy atom position relative to the
2265  // QM atom.
2266  bondVec *= bondVal ;
2267  }
2268 
2269  Position dummyPos = qmPos + bondVec;
2270 
2271  DebugM(1,"Creating dummy atom " << dummyPos << " ; QM ind: "
2272  << qmIt << " ; PC ind: " << pcIt << std::endl);
2273 
2274  dummyAtoms.push_back(dummyData(dummyPos, qmIt, qmMMBondIndx)) ;
2275 
2276  }
2277 
2278  DebugM(3, "Creating data for " << grpQMAtmVec.size() << " QM atoms "
2279  << dummyAtoms.size() << " dummy atoms " << grpPntChrgVec.size()
2280  << " real point charges and " << grpAppldChrgVec.size() - grpPntChrgVec.size()
2281  << " virtual point charges\n") ;
2282 
2283  int dataSize = grpQMAtmVec.size() + dummyAtoms.size() + grpAppldChrgVec.size();
2284  QMGrpCalcMsg *msg = new (dataSize, configLines.size(), 0) QMGrpCalcMsg;
2285  msg->timestep = timeStep;
2286  msg->grpIndx = grpIter;
2287  msg->peIter = peIter;
2288  msg->charge = grpChrg[grpIter];
2289  msg->multiplicity = grpMult[grpIter];
2290  msg->numQMAtoms = grpQMAtmVec.size();
2291  msg->numAllAtoms = grpQMAtmVec.size() + dummyAtoms.size();
2292  msg->numRealPntChrgs = grpPntChrgVec.size(); // The original set of classical atoms.
2293  msg->numAllPntChrgs = grpAppldChrgVec.size(); // The extended set with virtual point charges.
2294  msg->secProcOn = simParams->qmSecProcOn ;
2295  msg->constants = constants;
2296  msg->PMEOn = simParams->PMEOn ;
2297  if (msg->PMEOn)
2298  msg->PMEEwaldCoefficient = simParams->PMEEwaldCoefficient ;
2299  msg->switching = simParams->qmPCSwitchOn;
2300  msg->switchType = simParams->qmPCSwitchType;
2301  msg->cutoff = simParams->cutoff;
2302  msg->swdist = simParams->switchingDist;
2303  msg->pcScheme = simParams->qmPCScheme;
2304  msg->qmAtmChrgMode = simParams->qmChrgMode;
2305 
2306  strncpy(msg->baseDir, simParams->qmBaseDir, 256);
2307  strncpy(msg->execPath, simParams->qmExecPath, 256);
2308  if (msg->secProcOn)
2309  strncpy(msg->secProc, simParams->qmSecProc, 256);
2310 
2311  if (simParams->qmPrepProcOn && (timeStep == simParams->firstTimestep)) {
2312  msg->prepProcOn = true;
2313  strncpy(msg->prepProc, simParams->qmPrepProc, 256);
2314  } else
2315  msg->prepProcOn = false;
2316 
2317  QMAtomData *dataP = msg->data;
2318 
2319  for (int i=0; i<grpQMAtmVec.size(); i++) {
2320  dataP->position = grpQMAtmVec[i].position ;
2321  dataP->charge = grpQMAtmVec[i].charge ;
2322  dataP->id = grpQMAtmVec[i].id ;
2323  dataP->bountToIndx = -1;
2324  dataP->type = QMATOMTYPE_QM ;
2325  strncpy(dataP->element,elementArray[grpQMAtmVec[i].id],3);
2326  dataP++;
2327  }
2328 
2329  for (int i=0; i<dummyAtoms.size(); i++) {
2330  dataP->position = dummyAtoms[i].pos ;
2331  dataP->charge = 0 ;
2332  dataP->id = -1 ;
2333  dataP->bountToIndx = dummyAtoms[i].qmInd ;
2334  dataP->type = QMATOMTYPE_DUMMY ;
2335  strncpy(dataP->element,dummyElmArray[dummyAtoms[i].bondIndx],3);
2336  dataP++;
2337  }
2338 
2339  for (int i=0; i<grpAppldChrgVec.size(); i++) {
2340  dataP->position = grpAppldChrgVec[i].position ;
2341  dataP->charge = grpAppldChrgVec[i].charge ;
2342  // Point charges created to handle QM-MM bonds will have an id of -1.
2343  dataP->id = grpAppldChrgVec[i].id ;
2344  dataP->bountToIndx = -1;
2345  dataP->dist = grpAppldChrgVec[i].dist ;
2346  // If we are loading the classical atoms' charges
2347  // the point charge type is 1, unless it is from an
2348  // atom which is bound to a QM atom.
2349  if (i < grpPntChrgVec.size()) {
2350  if (grpAppldChrgVec[i].qmGrpID == -1) {
2351  dataP->type = QMPCTYPE_IGNORE ;
2352  }
2353  else if (grpAppldChrgVec[i].qmGrpID == -2) {
2354  dataP->type = QMPCTYPE_CLASSICAL ;
2355  dataP->bountToIndx = grpAppldChrgVec[i].homeIndx;
2356  }
2357  else {
2358  dataP->type = QMPCTYPE_CLASSICAL ;
2359  }
2360  }
2361  else {
2362  // Extra charges are created to handle QM-MM bonds (if they exist).
2363  dataP->type = QMPCTYPE_EXTRA ;
2364  dataP->bountToIndx = grpAppldChrgVec[i].id;
2365  }
2366  dataP++;
2367  }
2368 
2369  QMAtomData *qmP = msg->data ;
2370  QMAtomData *pcP = msg->data + msg->numAllAtoms ;
2371 
2372  // With this, every QM atom knows to which MM atom is is bound,
2373  // and vice versa. This will be usefull later on to prevent them from
2374  // feeling eachother's electrostatic charges AND to place the dummy
2375  // atom forces on the "real" atoms that form the bond.
2376  for( auto vecPtr = qmPCLocalIndxPairs.begin();
2377  vecPtr != qmPCLocalIndxPairs.end();
2378  vecPtr++ ) {
2379 
2380  int qmLocInd = (*vecPtr).first;
2381  int pcLocInd = (*vecPtr).second;
2382 
2383  qmP[qmLocInd].bountToIndx = pcLocInd ;
2384  pcP[pcLocInd].bountToIndx = qmLocInd ;
2385  }
2386 
2387 
2388  strcpy(msg->configLines, configLines.c_str());
2389 
2390  if (cSMDon)
2391  calcCSMD(grpIter, msg->numQMAtoms, qmP, cSMDForces) ;
2392 
2393  int targetPE = qmPEs[peIter] ;
2394 
2395  DebugM(4,"Sending QM group " << grpIter << " (ID " << grpID[grpIter]
2396  << ") to PE " << targetPE << std::endl);
2397 
2398  switch (simParams->qmFormat) {
2399  // Creates the command line that will be executed according to the
2400  // chosen QM software, as well as the input file with coordinates.
2401  case QMFormatORCA:
2402  QMProxy[targetPE].calcORCA(msg) ;
2403  break;
2404 
2405  case QMFormatMOPAC:
2406  QMProxy[targetPE].calcMOPAC(msg) ;
2407  break;
2408 
2409  case QMFormatUSR:
2410  QMProxy[targetPE].calcUSR(msg) ;
2411  break;
2412  }
2413 
2414  peIter++;
2415 
2416  if (peIter == qmPEs.size())
2417  peIter = 0;
2418  }
2419 
2420 }
Node::Object
static Node * Object()
Definition: Node.h:86
QMAtomData::bountToIndx
int bountToIndx
Definition: ComputeQM.C:299
QMGrpCalcMsg::PMEEwaldCoefficient
BigReal PMEEwaldCoefficient
Definition: ComputeQM.C:336
Molecule::get_qmMMNumTargs
const int * get_qmMMNumTargs()
Definition: Molecule.h:843
QMGrpCalcMsg::baseDir
char baseDir[256]
Definition: ComputeQM.C:338
QMRATIOTYPE
#define QMRATIOTYPE
Definition: ComputeQM.C:56
QMGrpCalcMsg::execPath
char execPath[256]
Definition: ComputeQM.C:338
QMGrpCalcMsg::data
QMAtomData * data
Definition: ComputeQM.C:339
QMGrpCalcMsg::multiplicity
Real multiplicity
Definition: ComputeQM.C:327
Molecule::get_qmGrpNumBonds
const int * get_qmGrpNumBonds()
Definition: Molecule.h:836
QMGrpCalcMsg::constants
BigReal constants
Definition: ComputeQM.C:328
QMAtomData::charge
float charge
Definition: ComputeQM.C:297
Vector
Definition: Vector.h:64
QMGrpCalcMsg::prepProcOn
bool prepProcOn
Definition: ComputeQM.C:330
ComputeQMAtom::charge
float charge
Definition: ComputeQM.C:67
QMATOMTYPE_DUMMY
#define QMATOMTYPE_DUMMY
Definition: ComputeQM.C:195
Real
float Real
Definition: common.h:109
COULOMB
#define COULOMB
Definition: common.h:46
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
QMGrpCalcMsg::numAllAtoms
int numAllAtoms
Definition: ComputeQM.C:324
Molecule::get_qmMMBondedIndx
const int *const * get_qmMMBondedIndx()
Definition: Molecule.h:839
QMPCTYPE_CLASSICAL
#define QMPCTYPE_CLASSICAL
Definition: ComputeQM.C:198
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
QMAtomData::dist
Real dist
Definition: ComputeQM.C:302
if
if(ComputeNonbondedUtil::goMethod==2)
Definition: ComputeNonbondedBase2.h:434
QMAtomData::position
Position position
Definition: ComputeQM.C:296
QMFormatMOPAC
#define QMFormatMOPAC
Definition: Molecule.h:130
iout
#define iout
Definition: InfoStream.h:51
Molecule::get_qmDummyElement
const char *const * get_qmDummyElement()
Definition: Molecule.h:840
QMGrpCalcMsg::switchType
int switchType
Definition: ComputeQM.C:333
QMFormatORCA
#define QMFormatORCA
Definition: Molecule.h:129
QMFormatUSR
#define QMFormatUSR
Definition: Molecule.h:131
Molecule::get_qmGrpBonds
const int *const * get_qmGrpBonds()
Definition: Molecule.h:838
QMGrpCalcMsg::switching
bool switching
Definition: ComputeQM.C:332
QMGrpCalcMsg::numQMAtoms
int numQMAtoms
Definition: ComputeQM.C:323
SortedArray::sort
void sort(void)
Definition: SortedArray.h:66
QMGrpCalcMsg::qmAtmChrgMode
int qmAtmChrgMode
Definition: ComputeQM.C:337
QMPCTYPE_IGNORE
#define QMPCTYPE_IGNORE
Definition: ComputeQM.C:197
QMGrpCalcMsg::numRealPntChrgs
int numRealPntChrgs
Definition: ComputeQM.C:325
Molecule::get_qmMMBond
const int *const * get_qmMMBond()
Definition: Molecule.h:837
QMGrpCalcMsg::prepProc
char prepProc[256]
Definition: ComputeQM.C:338
QMPCVec
std::vector< ComputeQMPntChrg > QMPCVec
Definition: ComputeQM.C:398
custom_ComputeQMAtom_Less
bool custom_ComputeQMAtom_Less(const ComputeQMAtom a, const ComputeQMAtom b)
Definition: ComputeQM.C:94
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:85
QMAtomData::element
char element[3]
Definition: ComputeQM.C:301
QMLENTYPE
#define QMLENTYPE
Definition: ComputeQM.C:55
Node::configList
ConfigList * configList
Definition: Node.h:179
Molecule::get_qmMMChargeTarget
const int *const * get_qmMMChargeTarget()
Definition: Molecule.h:842
Molecule::get_qmElements
const char *const * get_qmElements()
Definition: Molecule.h:827
QMATOMTYPE_QM
#define QMATOMTYPE_QM
Definition: ComputeQM.C:196
sort
BlockRadixSort::TempStorage sort
Definition: CudaTileListKernel.cu:693
simParams
#define simParams
Definition: Output.C:127
StringList::next
StringList * next
Definition: ConfigList.h:49
QMGrpCalcMsg::secProc
char secProc[256]
Definition: ComputeQM.C:338
StringList::data
char * data
Definition: ConfigList.h:48
ConfigList::find
StringList * find(const char *name) const
Definition: ConfigList.C:341
SortedArray::find
Elem * find(const Elem &elem)
QMForce::homeIndx
int homeIndx
Definition: ComputeQM.h:104
ResizeArray::size
int size(void) const
Definition: ResizeArray.h:127
iERROR
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:83
QMGrpCalcMsg::numAllPntChrgs
int numAllPntChrgs
Definition: ComputeQM.C:326
QMPCTYPE_EXTRA
#define QMPCTYPE_EXTRA
Definition: ComputeQM.C:199
QMGrpCalcMsg::configLines
char * configLines
Definition: ComputeQM.C:340
SortedArray::load
int load(const Elem &elem)
Definition: SortedArray.h:50
Molecule::get_qmDummyBondVal
const BigReal * get_qmDummyBondVal()
Definition: Molecule.h:835
QMGrpCalcMsg::secProcOn
bool secProcOn
Definition: ComputeQM.C:329
Vector::unit
Vector unit(void) const
Definition: Vector.h:182
BigReal
double BigReal
Definition: common.h:114
Molecule::get_qmGrpSizes
const int * get_qmGrpSizes()
Definition: Molecule.h:829
QMAtmVec
std::vector< ComputeQMAtom > QMAtmVec
Definition: ComputeQM.C:399
for
for(int i=0;i< n1;++i)
Definition: ComputeFullDirectBase.h:34
void ComputeQMMgr::setCompute ( ComputeQM *  c)
inline

Definition at line 406 of file ComputeQM.C.

406 { QMCompute = c; }
void ComputeQMMgr::storeQMRes ( QMGrpResMsg *  resMsg)

Definition at line 2422 of file ComputeQM.C.

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

2422  {
2423 
2424 // iout << iINFO << "Storing QM results for region " << resMsg->grpIndx
2425 // << " (ID: " << grpID[resMsg->grpIndx]
2426 // << ") with original energy: " << endi;
2427 // std::cout << std::fixed << std::setprecision(6) << resMsg->energyOrig << endi;
2428 // iout << " Kcal/mol\n" << endi;
2429 
2430 // if (resMsg->energyCorr != resMsg->energyOrig) {
2431 // iout << iINFO << "PME corrected energy: " << endi;
2432 // std::cout << std::fixed << std::setprecision(6) << resMsg->energyCorr << endi;
2433 // iout << " Kcal/mol\n" << endi;
2434 // }
2435 
2436  if ( (timeStep % simParams->qmEnergyOutFreq ) == 0 ) {
2437  iout << QMETITLE(timeStep);
2438  iout << FORMAT(grpID[resMsg->grpIndx]);
2439  iout << FORMAT(resMsg->energyOrig);
2440  if (resMsg->energyCorr != resMsg->energyOrig) iout << FORMAT(resMsg->energyCorr);
2441  iout << "\n\n" << endi;
2442  }
2443 
2444  totalEnergy += resMsg->energyCorr ;
2445 
2446  for ( int k=0; k<3; ++k )
2447  for ( int l=0; l<3; ++l )
2448  totVirial[k][l] += resMsg->virial[k][l];
2449 
2450  QMForce *fres = resMsg->force ;
2451  Real qmTotalCharge = 0;
2452 
2453  for (int i=0; i<resMsg->numForces; i++) {
2454 
2455  force[ fres[i].id ].force += fres[i].force;
2456 
2457  // Indicates the result is a QM atom, and we should get its updated charge.
2458  if (fres[i].replace == 1) {
2459  force[ fres[i].id ].charge = fres[i].charge;
2460  qmTotalCharge += fres[i].charge;
2461  }
2462  }
2463 
2464  if ((fabsf(roundf(qmTotalCharge) - qmTotalCharge) <= 0.001f) ) {
2465  qmTotalCharge = roundf(qmTotalCharge) ;
2466  }
2467 
2468  // In case we are calculating cSMD forces, apply them now.
2469  if (cSMDon) {
2470  for ( int i=0; i < cSMDindxLen[resMsg->grpIndx]; i++ ) {
2471  int cSMDit = cSMDindex[resMsg->grpIndx][i];
2472  force[ cSMDpairs[cSMDit][0] ].force += cSMDForces[cSMDit];
2473  }
2474  }
2475 
2476  DebugM(4,"QM total charge received is " << qmTotalCharge << std::endl);
2477 
2478  DebugM(4,"Current accumulated energy is " << totalEnergy << std::endl);
2479 
2480  numRecQMRes++;
2481 
2482  delete resMsg;
2483 }
QMGrpResMsg::force
QMForce * force
Definition: ComputeQM.C:181
Real
float Real
Definition: common.h:109
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
iout
#define iout
Definition: InfoStream.h:51
QMETITLE
static char * QMETITLE(int X)
Definition: ComputeQM.C:384
QMGrpResMsg::numForces
int numForces
Definition: ComputeQM.C:180
QMGrpResMsg::energyOrig
BigReal energyOrig
Definition: ComputeQM.C:177
QMForce::charge
float charge
Definition: ComputeQM.h:105
FORMAT
static char * FORMAT(BigReal X)
Definition: ComputeQM.C:367
QMForce::force
Force force
Definition: ComputeQM.h:103
simParams
#define simParams
Definition: Output.C:127
QMGrpResMsg::virial
BigReal virial[3][3]
Definition: ComputeQM.C:179
QMGrpResMsg::energyCorr
BigReal energyCorr
Definition: ComputeQM.C:178
for
for(int i=0;i< n1;++i)
Definition: ComputeFullDirectBase.h:34
The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.5