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

Constructor & Destructor Documentation

◆ ComputeQMMgr()

ComputeQMMgr::ComputeQMMgr ( )

Definition at line 549 of file ComputeQM.C.

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

◆ ~ComputeQMMgr()

ComputeQMMgr::~ComputeQMMgr ( )

Definition at line 558 of file ComputeQM.C.

References close_dcd_write().

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

Member Function Documentation

◆ calcMOPAC()

void ComputeQMMgr::calcMOPAC ( QMGrpCalcMsg *  msg)

Definition at line 2802 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, 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, Vector::y, and Vector::z.

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

◆ calcORCA()

void ComputeQMMgr::calcORCA ( QMGrpCalcMsg *  msg)

Definition at line 3650 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, 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, Vector::y, and Vector::z.

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

◆ calcUSR()

void ComputeQMMgr::calcUSR ( QMGrpCalcMsg *  msg)

Definition at line 4512 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, 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, Vector::y, and Vector::z.

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

◆ get_subsArray()

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

Definition at line 429 of file ComputeQM.C.

Referenced by lssSubs().

429 { return subsArray; } ;

◆ procQMRes()

void ComputeQMMgr::procQMRes ( )

Definition at line 2486 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< Elem >::insert(), iout, QMForceMsg::lssDat, QMCoordMsg::numAtoms, QMPntChrgMsg::numAtoms, QMForceMsg::numForces, QMForceMsg::numLssDat, QMForceMsg::PMEOn, SimParameters::PMEOn, ComputeQMAtom::position, SimParameters::qmOutFreq, SimParameters::qmPosOutFreq, ResizeArray< Elem >::size(), SortedArray< Elem >::sort(), QMCoordMsg::sourceNode, QMPntChrgMsg::sourceNode, QMForceMsg::virial, write_dcdstep(), Vector::x, Vector::y, and Vector::z.

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

◆ recvForce()

void ComputeQMMgr::recvForce ( QMForceMsg *  fmsg)

Definition at line 2663 of file ComputeQM.C.

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

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

◆ recvFullQM()

void ComputeQMMgr::recvFullQM ( QMCoordMsg *  qmFullMsg)

Definition at line 1282 of file ComputeQM.C.

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

1282  {
1283 
1284  if (subsArray.size() > 0)
1285  subsArray.clear();
1286 
1287  QMCompute->processFullQM(qmFullMsg);
1288 }
ResizeArray::size
int size(void) const
Definition: ResizeArray.h:131
ResizeArray::clear
void clear()
Definition: ResizeArray.h:91
ComputeQM::processFullQM
void processFullQM(QMCoordMsg *)
Definition: ComputeQM.C:1293

◆ recvPartQM()

void ComputeQMMgr::recvPartQM ( QMCoordMsg *  msg)

Definition at line 819 of file ComputeQM.C.

References ResizeArray< Elem >::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< Elem >::size(), QMCoordMsg::sourceNode, QMPntChrgMsg::sourceNode, TIMEFACTOR, QMCoordMsg::timestep, ComputeQMPntChrg::vdwType, and write_dcdheader().

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

◆ recvPntChrg()

void ComputeQMMgr::recvPntChrg ( QMPntChrgMsg *  msg)

Definition at line 1911 of file ComputeQM.C.

References QMGrpCalcMsg::baseDir, QMAtomData::bountToIndx, ComputeQMAtom::charge, QMAtomData::charge, QMGrpCalcMsg::charge, ResizeArray< Elem >::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< Elem >::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< Elem >::load(), QMGrpCalcMsg::multiplicity, NAMD_die(), StringList::next, SimParameters::nonbondedScaling, QMGrpCalcMsg::numAllAtoms, QMGrpCalcMsg::numAllPntChrgs, QMPntChrgMsg::numAtoms, QMGrpCalcMsg::numQMAtoms, QMGrpCalcMsg::numRealPntChrgs, Node::Object(), QMGrpCalcMsg::pcScheme, QMGrpCalcMsg::peIter, 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< Elem >::size(), SortedArray< Elem >::sort(), QMGrpCalcMsg::swdist, QMGrpCalcMsg::switching, SimParameters::switchingDist, QMGrpCalcMsg::switchType, QMGrpCalcMsg::timestep, QMAtomData::type, and Vector::unit().

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

◆ setCompute()

void ComputeQMMgr::setCompute ( ComputeQM *  c)
inline

Definition at line 407 of file ComputeQM.C.

407 { QMCompute = c; }

◆ storeQMRes()

void ComputeQMMgr::storeQMRes ( QMGrpResMsg *  resMsg)

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

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