C++ Lattice::b_p方法代码示例

void ComputeExtMgr::recvCoord(ExtCoordMsg *msg) {
  if ( ! numSources ) {
    numSources = (PatchMap::Object())->numNodesWithPatches();
    coordMsgs = new ExtCoordMsg*[numSources];
    for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
    numArrived = 0;
    numAtoms = Node::Object()->molecule->numAtoms;
    coord = new ComputeExtAtom[numAtoms];
    force = new ExtForce[numAtoms];
  }

  int i;
  for ( i=0; i < msg->numAtoms; ++i ) {
    coord[msg->coord[i].id] = msg->coord[i];
  }

  coordMsgs[numArrived] = msg;
  ++numArrived;

  if ( numArrived < numSources ) return;
  numArrived = 0;

  // ALL DATA ARRIVED --- CALCULATE FORCES
  Lattice lattice = msg->lattice;
  SimParameters *simParams = Node::Object()->simParameters;
  FILE *file;
  int iret;

  // write coordinates to file
  //iout << "writing to file " << simParams->extCoordFilename << "\n" << endi;
  file = fopen(simParams->extCoordFilename,"w");
  if ( ! file ) { NAMD_die(strerror(errno)); }
  for ( i=0; i<numAtoms; ++i ) {
    int id = coord[i].id + 1;
    double charge = coord[i].charge;
    double x = coord[i].position.x;
    double y = coord[i].position.y;
    double z = coord[i].position.z;
    iret = fprintf(file,"%d %f %f %f %f\n",id,charge,x,y,z);
    if ( iret < 0 ) { NAMD_die(strerror(errno)); }
  }
  // write periodic cell lattice (0 0 0 if non-periodic)
  Vector a = lattice.a();  if ( ! lattice.a_p() ) a = Vector(0,0,0);
  Vector b = lattice.b();  if ( ! lattice.b_p() ) b = Vector(0,0,0);
  Vector c = lattice.c();  if ( ! lattice.c_p() ) c = Vector(0,0,0);
  iret = fprintf(file,"%f %f %f\n%f %f %f\n%f %f %f\n",
                 a.x, a.y, a.z, b.x, b.y, b.z, c.x, c.y, c.z);
  if ( iret < 0 ) { NAMD_die(strerror(errno)); }
  fclose(file);

  // run user-specified command
  //iout << "running command " << simParams->extForcesCommand << "\n" << endi;
  iret = system(simParams->extForcesCommand);
  if ( iret == -1 ) { NAMD_die(strerror(errno)); }
  if ( iret ) { NAMD_die("Error running command for external forces."); }

  // remove coordinate file
  iret = remove(simParams->extCoordFilename);
  if ( iret ) { NAMD_die(strerror(errno)); }

  // read forces from file (overwrite positions)
  //iout << "reading from file " << simParams->extForceFilename << "\n" << endi;
  file = fopen(simParams->extForceFilename,"r");
  if ( ! file ) { NAMD_die(strerror(errno)); }
  for ( i=0; i<numAtoms; ++i ) {
    int id, replace;
    double x, y, z;
    iret = fscanf(file,"%d %d %lf %lf %lf\n", &id, &replace, &x, &y, &z);
    if ( iret != 5 ) { NAMD_die("Error reading external forces file."); }
    if ( id != i + 1 ) { NAMD_die("Atom ID error in external forces file."); }
    force[i].force.x = x; force[i].force.y = y; force[i].force.z = z;
    force[i].replace = replace;
  }
  // read energy and virial if they are present
  // virial used by NAMD is -'ve of normal convention, so reverse it!
  // virial[i][j] in file should be sum of -1 * f_i * r_j
  double energy;
  double virial[3][3];
  iret = fscanf(file,"%lf\n", &energy);
  if ( iret != 1 ) {
    energy = 0;
    for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] = 0;
  } else {
    iret = fscanf(file,"%lf %lf %lf\n%lf %lf %lf\n%lf %lf %lf\n",
      &virial[0][0], &virial[0][1], &virial[0][2],
      &virial[1][0], &virial[1][1], &virial[1][2],
      &virial[2][0], &virial[2][1], &virial[2][2]);
    if ( iret != 9 ) {
      for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] = 0;
    } else {
      // virial used by NAMD is -'ve of normal convention, so reverse it!
      for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] *= -1.0;
    }
  }
  fclose(file);

  // remove force file
  iret = remove(simParams->extForceFilename);
  if ( iret ) { NAMD_die(strerror(errno)); }

//.........这里部分代码省略.........

void ComputeMsmSerialMgr::recvCoord(MsmSerialCoordMsg *msg) {
  if ( ! numSources ) {
    numSources = (PatchMap::Object())->numNodesWithPatches();
    coordMsgs = new MsmSerialCoordMsg*[numSources];
    for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
    numArrived = 0;
    numAtoms = Node::Object()->molecule->numAtoms;
    coord = new ComputeMsmSerialAtom[numAtoms];
    force = new MsmSerialForce[numAtoms];
  }

  int i;
  for ( i=0; i < msg->numAtoms; ++i ) {
    coord[msg->coord[i].id] = msg->coord[i];
  }

  coordMsgs[numArrived] = msg;
  ++numArrived;

  if ( numArrived < numSources ) return;
  numArrived = 0;

  // ALL DATA ARRIVED --- CALCULATE FORCES
  Lattice lattice = msg->lattice;
  SimParameters *simParams = Node::Object()->simParameters;

  double energy = 0;
  double virial[3][3];

  int rc = 0;  // return code

  if ( ! msmsolver ) {
    //
    // setup MSM solver
    //
    msmsolver = NL_msm_create();
    if ( ! msmsolver ) NAMD_die("unable to create MSM solver");
    double dielectric = simParams->dielectric;
    double cutoff = simParams->cutoff;
    double gridspacing = simParams->MSMGridSpacing;
    double padding = simParams->MSMPadding;
    int approx = simParams->MSMApprox;
    int split = simParams->MSMSplit;
    int nlevels = simParams->MSMLevels;
    int msmflags = 0;
    msmflags |= (lattice.a_p() ? NL_MSM_PERIODIC_VEC1 : 0);
    msmflags |= (lattice.b_p() ? NL_MSM_PERIODIC_VEC2 : 0);
    msmflags |= (lattice.c_p() ? NL_MSM_PERIODIC_VEC3 : 0);
    msmflags |= NL_MSM_COMPUTE_LONG_RANGE;  // compute only long-range part
    //msmflags |= NL_MSM_COMPUTE_ALL;
    //printf("msmflags = %x\n", msmflags);
    rc = NL_msm_configure(msmsolver, gridspacing, approx, split, nlevels);
    if (rc) NAMD_die("unable to configure MSM solver");
    Vector u=lattice.a(), v=lattice.b(), w=lattice.c(), c=lattice.origin();
    Vector ru=lattice.a_r(), rv=lattice.b_r(), rw=lattice.c_r();
    if ((msmflags & NL_MSM_PERIODIC_ALL) != NL_MSM_PERIODIC_ALL) {
      // called only if there is some non-periodic boundary
      int isperiodic = (msmflags & NL_MSM_PERIODIC_ALL);
      //printf("calling rescale\n");
      rescale_nonperiodic_cell(u, v, w, c, ru, rv, rw,
          isperiodic, numAtoms, coord, padding, gridspacing);
    }
    double vec1[3], vec2[3], vec3[3], center[3];
    vec1[0] = u.x;
    vec1[1] = u.y;
    vec1[2] = u.z;
    vec2[0] = v.x;
    vec2[1] = v.y;
    vec2[2] = v.z;
    vec3[0] = w.x;
    vec3[1] = w.y;
    vec3[2] = w.z;
    center[0] = c.x;
    center[1] = c.y;
    center[2] = c.z;
#if 0
    printf("dielectric = %g\n", dielectric);
    printf("vec1 = %g %g %g\n", vec1[0], vec1[1], vec1[2]);
    printf("vec2 = %g %g %g\n", vec2[0], vec2[1], vec2[2]);
    printf("vec3 = %g %g %g\n", vec3[0], vec3[1], vec3[2]);
    printf("center = %g %g %g\n", center[0], center[1], center[2]);
    printf("cutoff = %g\n", cutoff);
    printf("numatoms = %d\n", numAtoms);
#endif
    rc = NL_msm_setup(msmsolver, cutoff, vec1, vec2, vec3, center, msmflags);
    if (rc) NAMD_die("unable to set up MSM solver");
    msmcoord = new double[4*numAtoms];
    msmforce = new double[3*numAtoms];
    if (msmcoord==0 || msmforce==0) NAMD_die("can't allocate MSM atom buffers");
    // scale charges - these won't change
    double celec = sqrt(COULOMB / simParams->dielectric);
    for (i = 0;  i < numAtoms;  i++) {
      msmcoord[4*i+3] = celec * coord[i].charge;
    }
  }

  // evaluate long-range MSM forces
  for (i = 0;  i < numAtoms;  i++) {
    msmcoord[4*i  ] = coord[i].position.x;
    msmcoord[4*i+1] = coord[i].position.y;
//.........这里部分代码省略.........

int ScriptTcl::Tcl_coorfile(ClientData clientData,
	Tcl_Interp *interp, int argc, char *argv[]) {
  ScriptTcl *script = (ScriptTcl *)clientData;
  script->initcheck();
  if (argc == 4 && !strcmp(argv[1], "open")) {
    if (strcmp(argv[2], "dcd")) {
      NAMD_die("Sorry, coorfile presently supports only DCD files");
    }
    filehandle = dcdplugin->open_file_read(argv[3], "dcd", &numatoms);
    if (!filehandle) {
      Tcl_AppendResult(interp, "coorfile: Error opening file ", argv[3], NULL);
      return TCL_ERROR;
    }
    if (numatoms != Node::Object()->pdb->num_atoms()) {
      Tcl_AppendResult(interp, "Coordinate file ", argv[3],
        "\ncontains the wrong number of atoms.", NULL);
      return TCL_ERROR;
    }
    coords = new float[3*numatoms];
    vcoords = new Vector[3*numatoms];
    iout << iINFO << "Coordinate file " << argv[3] << " opened for reading.\n"
         << endi;
  } else if (argc == 2 && !strcmp(argv[1], "read")) {
    if (filehandle == NULL) {
      Tcl_AppendResult(interp, "coorfile read: Error, no file open for reading",
	NULL);
      return TCL_ERROR;
    }
    molfile_timestep_t ts;
    ts.coords = coords;
    int rc = dcdplugin->read_next_timestep(filehandle, numatoms, &ts);
    if (rc) {  // EOF
      Tcl_SetObjResult(interp, Tcl_NewIntObj(-1));
      return TCL_OK;
    }
    iout << iINFO << "Reading timestep from file.\n" << endi;
    Lattice lattice;
    if (get_lattice_from_ts(&lattice, &ts)) {
      iout << iINFO << "Updating unit cell from timestep.\n" << endi;
      if ( lattice.a_p() && ! script->state->lattice.a_p() ||
           lattice.b_p() && ! script->state->lattice.b_p() ||
           lattice.c_p() && ! script->state->lattice.c_p() ) {
        iout << iWARN << "Cell basis vectors should be specified before reading trajectory.\n" << endi;
      }
      // update Controller's lattice, but don't change the origin!
      Vector a(0.);  if ( script->state->lattice.a_p() ) a = lattice.a();
      Vector b(0.);  if ( script->state->lattice.b_p() ) b = lattice.b();
      Vector c(0.);  if ( script->state->lattice.c_p() ) c = lattice.c();
      script->state->lattice.set(a,b,c);
      SetLatticeMsg *msg = new SetLatticeMsg;
      msg->lattice = script->state->lattice;
      (CProxy_PatchMgr(CkpvAccess(BOCclass_group).patchMgr)).setLattice(msg);
      script->barrier();
    }
    for (int i=0; i<numatoms; i++) {
      vcoords[i].x = coords[3*i+0];
      vcoords[i].y = coords[3*i+1];
      vcoords[i].z = coords[3*i+2];
    }
    Node::Object()->pdb->set_all_positions(vcoords);
    script->reinitAtoms();
    Tcl_SetObjResult(interp, Tcl_NewIntObj(0));
  } else if (argc == 2 && !strcmp(argv[1], "close")) {
    if (!filehandle) {
      Tcl_AppendResult(interp, "coorfile close: No file opened for reading!",
        NULL);
      return TCL_OK;
    }
    iout << iINFO << "Closing coordinate file.\n" << endi;
    dcdplugin->close_file_read(filehandle);
    filehandle = NULL;
    delete [] coords;
    delete [] vcoords;

  } else if (argc ==2 && !strcmp(argv[1], "skip")) {
    if (filehandle == NULL) {
      Tcl_AppendResult(interp, "coorfile skip: Error, no file open for reading",
	NULL);
      return TCL_ERROR;
    }
    int rc = dcdplugin->read_next_timestep(filehandle, numatoms, NULL);
    if (rc) {  // EOF
      Tcl_SetObjResult(interp, Tcl_NewIntObj(-1));
      return TCL_OK;
    }
    Tcl_SetObjResult(interp, Tcl_NewIntObj(0));

  } else {
    NAMD_die("Unknown option passed to coorfile");
  }
  return TCL_OK;
}