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

double PmeKSpace::compute_energy(float *q_arr, const Lattice &lattice, double ewald, double *virial) {
  double energy = 0.0;
  double v0 = 0.;
  double v1 = 0.;
  double v2 = 0.;
  double v3 = 0.;
  double v4 = 0.;
  double v5 = 0.;

  int n;
  int k1, k2, k3, ind;
  int K1, K2, K3;

  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3;

  i_pi_volume = 1.0/(M_PI * lattice.volume());
  piob = M_PI/ewald;
  piob *= piob;

  if ( lattice.orthogonal() ) {
  // if ( 0 ) { // JCP FOR TESTING
    //This branch is the usual call path.
#if     USE_CKLOOP
    int useCkLoop = Node::Object()->simParameters->useCkLoop;
    if(useCkLoop>=CKLOOP_CTRL_PME_KSPACE){
        return compute_energy_orthogonal_helper(q_arr, lattice, ewald, virial);
    }
#endif    
    double recipx = lattice.a_r().x;
    double recipy = lattice.b_r().y;
    double recipz = lattice.c_r().z;
        
    init_exp(exp1, K1, 0, K1, recipx);
    init_exp(exp2, K2, k2_start, k2_end, recipy);
    init_exp(exp3, K3, k3_start, k3_end, recipz);

    ind = 0;
    for ( k1=0; k1<K1; ++k1 ) {
      double m1, m11, b1, xp1;
      b1 = bm1[k1];
      int k1_s = k1<=K1/2 ? k1 : k1-K1;
      m1 = k1_s*recipx;
      m11 = m1*m1;
      xp1 = i_pi_volume*exp1[abs(k1_s)];
      for ( k2=k2_start; k2<k2_end; ++k2 ) {
        double m2, m22, b1b2, xp2;
        b1b2 = b1*bm2[k2];
        int k2_s = k2<=K2/2 ? k2 : k2-K2;
        m2 = k2_s*recipy;
        m22 = m2*m2;
        xp2 = exp2[abs(k2_s)]*xp1;
        k3 = k3_start;
        if ( k1==0 && k2==0 && k3==0 ) {
          q_arr[ind++] = 0.0;
          q_arr[ind++] = 0.0;
          ++k3;
        }
        for ( ; k3<k3_end; ++k3 ) {
          double m3, m33, xp3, msq, imsq, vir, fac;
          double theta3, theta, q2, qr, qc, C;
          theta3 = bm3[k3] *b1b2;
          m3 = k3*recipz;
          m33 = m3*m3;
          xp3 = exp3[k3];
          qr = q_arr[ind]; qc=q_arr[ind+1];
          q2 = 2*(qr*qr + qc*qc)*theta3;
          if ( (k3 == 0) || ( k3 == K3/2 && ! (K3 & 1) ) ) q2 *= 0.5;
          msq = m11 + m22 + m33;
          imsq = 1.0/msq;
          C = xp2*xp3*imsq;
          theta = theta3*C;
          q_arr[ind] *= theta;
          q_arr[ind+1] *= theta;
          vir = -2*(piob+imsq);
          fac = q2*C;
          energy += fac;
          v0 += fac*(1.0+vir*m11);
          v1 += fac*vir*m1*m2;
          v2 += fac*vir*m1*m3;
          v3 += fac*(1.0+vir*m22);
          v4 += fac*vir*m2*m3;
          v5 += fac*(1.0+vir*m33);
          ind += 2;
        }
      }
    }
    
  } else if ( cross(lattice.a(),lattice.b()).unit() == lattice.c().unit() ) {
  // } else if ( 0 ) { // JCP FOR TESTING
    Vector recip1 = lattice.a_r();
    Vector recip2 = lattice.b_r();
    Vector recip3 = lattice.c_r();
    double recip3_x = recip3.x;
    double recip3_y = recip3.y;
    double recip3_z = recip3.z;
    init_exp(exp3, K3, k3_start, k3_end, recip3.length());

    ind = 0;
    for ( k1=0; k1<K1; ++k1 ) {
      double b1; Vector m1;
//.........这里部分代码省略.........

void ComputeDPMEMaster::recvData(ComputeDPMEDataMsg *msg)
{ 
  DebugM(4,"ComputeDPMEMaster::recvData() " << msg->numParticles
	<< " particles from node " << msg->node << "\n");

  {
    homeNode.add(msg->node);
    Pme2Particle *data_ptr = localData + numLocalAtoms;
    for ( int j = 0; j < msg->numParticles; ++j, ++data_ptr ) {
      *data_ptr = msg->particles[j];
    }
    numLocalAtoms += msg->numParticles;
    endForNode.add(numLocalAtoms);
    delete msg;
  }

  if ( homeNode.size() < numWorkingPes ) return;  // messages outstanding

  DebugM(4,"ComputeDPMEMaster::recvData() running serial code.\n");

  // single processor version

  Lattice lattice = host->getFlags()->lattice;
  SimParameters * simParams = Node::Object()->simParameters;
  int i;

  AtomInfo atom_info;
  atom_info.numatoms = numLocalAtoms;
  atom_info.atompnt = 0;  // not used
  atom_info.freepnt = 0;  // not used
  atom_info.nlocal = numLocalAtoms;
  atom_info.nother = 0;

  if ( ! lattice.orthogonal() ) {
    NAMD_die("DPME only supports orthogonal PBC's.");
  }

  BoxInfo box_info;
  box_info.box.x = lattice.a().x;
  box_info.box.y = lattice.b().y;
  box_info.box.z = lattice.c().z;
  box_info.box.origin = 0.;  // why only one number?
  box_info.prd.x = box_info.box.x;
  box_info.prd.y = box_info.box.y;
  box_info.prd.z = box_info.box.z;
  box_info.prd2.x = 0.5 * box_info.box.x;
  box_info.prd2.y = 0.5 * box_info.box.y;
  box_info.prd2.z = 0.5 * box_info.box.z;
  box_info.cutoff = simParams->cutoff;
  box_info.rs = simParams->cutoff;
  box_info.mc2.x = 2. * ( box_info.prd.x - box_info.cutoff );
  box_info.mc2.y = 2. * ( box_info.prd.y - box_info.cutoff );
  box_info.mc2.z = 2. * ( box_info.prd.z - box_info.cutoff );
  box_info.ewaldcof = ComputeNonbondedUtil::ewaldcof;
  box_info.dtol = simParams->PMETolerance;
  for (i = 0; i <= 8; i++) {
    box_info.recip[i] = 0.; /* assume orthogonal box */
  }
  box_info.recip[0] = 1.0/box_info.box.x;
  box_info.recip[4] = 1.0/box_info.box.y;
  box_info.recip[8] = 1.0/box_info.box.z;

  GridInfo grid_info;
  grid_info.order = simParams->PMEInterpOrder;
  grid_info.nfftgrd.x = simParams->PMEGridSizeX;
  grid_info.nfftgrd.y = simParams->PMEGridSizeY;
  grid_info.nfftgrd.z = simParams->PMEGridSizeZ;
  grid_info.nfft = 0;
  grid_info.volume = lattice.volume();

  PeInfo pe_info;  // hopefully this isn't used for anything
  pe_info.myproc.node = 0;
  pe_info.myproc.nprocs = 1;
  pe_info.myproc.ncube = 0;
  pe_info.inst_node[0] = 0;
  pe_info.igrid = 0;

  PmeVector *localResults;
  double recip_vir[6];
  int time_count = 0;
  int tsteps = 1;
  double mytime = 0.;

  // perform calculations
  BigReal electEnergy = 0;

  // calculate self energy
  Pme2Particle *data_ptr = localData;
  for(i=0; i<numLocalAtoms; ++i)
  {
    electEnergy += data_ptr->cg * data_ptr->cg;
    ++data_ptr;
  }
  electEnergy *= -1. * box_info.ewaldcof / SQRT_PI;

  DebugM(4,"Ewald self energy: " << electEnergy << "\n");

  DebugM(4,"Calling dpme_eval_recip().\n");

  double pme_start_time = 0;
//.........这里部分代码省略.........

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;
}