C++ pr_rvecs函数代码示例

void list_xtc(const char *fn)
{
    t_fileio  *xd;
    int        indent;
    char       buf[256];
    rvec      *x;
    matrix     box;
    int        nframe, natoms, step;
    real       prec, time;
    gmx_bool   bOK;

    xd = open_xtc(fn, "r");
    read_first_xtc(xd, &natoms, &step, &time, box, &x, &prec, &bOK);

    nframe = 0;
    do
    {
        sprintf(buf, "%s frame %d", fn, nframe);
        indent = 0;
        indent = pr_title(stdout, indent, buf);
        pr_indent(stdout, indent);
        fprintf(stdout, "natoms=%10d  step=%10d  time=%12.7e  prec=%10g\n",
                natoms, step, time, prec);
        pr_rvecs(stdout, indent, "box", box, DIM);
        pr_rvecs(stdout, indent, "x", x, natoms);
        nframe++;
    }
    while (read_next_xtc(xd, natoms, &step, &time, box, x, &prec, &bOK));
    if (!bOK)
    {
        fprintf(stderr, "\nWARNING: Incomplete frame at time %g\n", time);
    }
    sfree(x);
    close_xtc(xd);
}

/* Check whether the box is unchanged */
static void check_box_unchanged(matrix f_box, matrix box, char fn[], warninp_t wi)
{
    int      i, ii;
    gmx_bool bSame = TRUE;
    char     warn_buf[STRLEN];


    for (i = 0; i < DIM; i++)
    {
        for (ii = 0; ii < DIM; ii++)
        {
            if (f_box[i][ii] != box[i][ii])
            {
                bSame = FALSE;
            }
        }
    }
    if (!bSame)
    {
        sprintf(warn_buf, "%s Box size in reference file %s differs from actual box size!",
                RotStr, fn);
        warning(wi, warn_buf);
        pr_rvecs(stderr, 0, "Your box is:", box, 3);
        pr_rvecs(stderr, 0, "Box in file:", f_box, 3);
    }
}

static void sum_forces(int start,int end,rvec f[],rvec flr[])
{
  int i;
  
  if (gmx_debug_at) {
    pr_rvecs(debug,0,"fsr",f+start,end-start);
    pr_rvecs(debug,0,"flr",flr+start,end-start);
  }
  for(i=start; (i<end); i++)
    rvec_inc(f[i],flr[i]);
}

static void list_trn(char *fn)
{
  int         fpread,fpwrite,nframe,indent;
  char        buf[256];
  rvec        *x,*v,*f;
  matrix      box;
  t_trnheader trn;
  bool        bOK;

  fpread  = open_trn(fn,"r"); 
  fpwrite = open_tpx(NULL,"w");
  gmx_fio_setdebug(fpwrite,TRUE);
  
  nframe = 0;
  while (fread_trnheader(fpread,&trn,&bOK)) {
    snew(x,trn.natoms);
    snew(v,trn.natoms);
    snew(f,trn.natoms);
    if (fread_htrn(fpread,&trn,
		   trn.box_size ? box : NULL,
		   trn.x_size   ? x : NULL,
		   trn.v_size   ? v : NULL,
		   trn.f_size   ? f : NULL)) {
      sprintf(buf,"%s frame %d",fn,nframe);
      indent=0;
      indent=pr_title(stdout,indent,buf);
      pr_indent(stdout,indent);
      fprintf(stdout,"natoms=%10d  step=%10d  time=%12.7e  lambda=%10g\n",
	      trn.natoms,trn.step,trn.t,trn.lambda);
      if (trn.box_size)
	pr_rvecs(stdout,indent,"box",box,DIM);
      if (trn.x_size)
	pr_rvecs(stdout,indent,"x",x,trn.natoms);
      if (trn.v_size)
	pr_rvecs(stdout,indent,"v",v,trn.natoms);
      if (trn.f_size)
	pr_rvecs(stdout,indent,"f",f,trn.natoms);
    } 
    else
      fprintf(stderr,"\nWARNING: Incomplete frame: nr %d, t=%g\n",
	      nframe,trn.t);
    
    sfree(x);
    sfree(v);
    sfree(f);
    nframe++;
  }
  if (!bOK)
    fprintf(stderr,"\nWARNING: Incomplete frame header: nr %d, t=%g\n",
	    nframe,trn.t);
  close_tpx(fpwrite);
  close_trn(fpread);
}

static void dump_it_all(FILE *fp,char *title,
                        int natoms,rvec x[],rvec xp[],rvec v[],rvec f[])
{
#ifdef DEBUG
  if (fp) {
    fprintf(fp,"%s\n",title);
    pr_rvecs(fp,0,"x",x,natoms);
    pr_rvecs(fp,0,"xp",xp,natoms);
    pr_rvecs(fp,0,"v",v,natoms);
    pr_rvecs(fp,0,"f",f,natoms);
  }
#endif
}

static int correct_box_elem(FILE *fplog, int step, tensor box, int v, int d)
{
    int shift, maxshift = 10;

    shift = 0;

    /* correct elem d of vector v with vector d */
    while (box[v][d] > BOX_MARGIN_CORRECT*0.5*box[d][d])
    {
        if (fplog)
        {
            fprintf(fplog, "Step %d: correcting invalid box:\n", step);
            pr_rvecs(fplog, 0, "old box", box, DIM);
        }
        rvec_dec(box[v], box[d]);
        shift--;
        if (fplog)
        {
            pr_rvecs(fplog, 0, "new box", box, DIM);
        }
        if (shift <= -maxshift)
        {
            gmx_fatal(FARGS,
                      "Box was shifted at least %d times. Please see log-file.",
                      maxshift);
        }
    }
    while (box[v][d] < -BOX_MARGIN_CORRECT*0.5*box[d][d])
    {
        if (fplog)
        {
            fprintf(fplog, "Step %d: correcting invalid box:\n", step);
            pr_rvecs(fplog, 0, "old box", box, DIM);
        }
        rvec_inc(box[v], box[d]);
        shift++;
        if (fplog)
        {
            pr_rvecs(fplog, 0, "new box", box, DIM);
        }
        if (shift >= maxshift)
        {
            gmx_fatal(FARGS,
                      "Box was shifted at least %d times. Please see log-file.",
                      maxshift);
        }
    }

    return shift;
}

void calc_force(int natom,rvec f[],rvec fff[])
{
    int  i,j,m;
    int  jindex[] = { 0, 5, 10};
    rvec dx,df;
    real msf1,msf2;

    for(j=0; (j<2); j++) {
        clear_rvec(fff[j]);
        for(i=jindex[j]; (i<jindex[j+1]); i++) {
            for(m=0; (m<DIM); m++) {
                fff[j][m] += f[i][m];
            }
        }
    }

    msf1 = iprod(fff[0],fff[0]);
    msf2 = iprod(fff[1],fff[1]);
    if (debug) {
        pr_rvecs(debug,0,"force",f,natom);

        fprintf(debug,"FMOL:  %10.3f  %10.3f  %10.3f  %10.3f  %10.3f  %10.3f\n",
                fff[0][XX],fff[0][YY],fff[0][ZZ],fff[1][XX],fff[1][YY],fff[1][ZZ]);
        fprintf(debug,"RMSF:  %10.3e  %10.3e\n",msf1,msf2);
    }
}

static void calc_virial(FILE *fplog,int start,int homenr,rvec x[],rvec f[],
			tensor vir_part,t_graph *graph,matrix box,
			t_nrnb *nrnb,const t_forcerec *fr,int ePBC)
{
  int i,j;
  tensor virtest;

  /* The short-range virial from surrounding boxes */
  clear_mat(vir_part);
  calc_vir(fplog,SHIFTS,fr->shift_vec,fr->fshift,vir_part,ePBC==epbcSCREW,box);
  inc_nrnb(nrnb,eNR_VIRIAL,SHIFTS);
  
  /* Calculate partial virial, for local atoms only, based on short range. 
   * Total virial is computed in global_stat, called from do_md 
   */
  f_calc_vir(fplog,start,start+homenr,x,f,vir_part,graph,box);
  inc_nrnb(nrnb,eNR_VIRIAL,homenr);

  /* Add position restraint contribution */
  for(i=0; i<DIM; i++) {
    vir_part[i][i] += fr->vir_diag_posres[i];
  }

  /* Add wall contribution */
  for(i=0; i<DIM; i++) {
    vir_part[i][ZZ] += fr->vir_wall_z[i];
  }

  if (debug)
    pr_rvecs(debug,0,"vir_part",vir_part,DIM);
}

static void pr_matrix(FILE *fp,int indent,const char *title,rvec *m,
		      gmx_bool bMDPformat)
{
  if (bMDPformat)
    fprintf(fp,"%-10s    = %g %g %g %g %g %g\n",title,
	    m[XX][XX],m[YY][YY],m[ZZ][ZZ],m[XX][YY],m[XX][ZZ],m[YY][ZZ]);
  else
    pr_rvecs(fp,indent,title,m,DIM);
}

void init_parallel(FILE *log,char *tpxfile,t_commrec *cr,
		   t_inputrec *inputrec,gmx_mtop_t *mtop,
		   t_state *state,
		   int list)
{
  int  step;
  real t;
  char buf[256];
  
  if (MASTER(cr)) {
    init_inputrec(inputrec);
    read_tpx_state(tpxfile,&step,&t,inputrec,state,NULL,mtop);
    /* When we will be doing domain decomposition with separate PME nodes
     * the rng entries will be too large, we correct for this later.
     */
    set_state_entries(state,inputrec,cr->nnodes);
  }
  bcast_ir_mtop(cr,inputrec,mtop);

  if (inputrec->eI == eiBD || EI_SD(inputrec->eI)) {
    /* Make sure the random seeds are different on each node */
    inputrec->ld_seed += cr->nodeid;
  }
  
  /* Printing */
  if (list!=0 && log!=NULL) 
  {
	  if (list&LIST_INPUTREC)
		  pr_inputrec(log,0,"parameters of the run",inputrec,FALSE);
	  if (list&LIST_X)
		  pr_rvecs(log,0,"box",state->box,DIM);
	  if (list&LIST_X)
		  pr_rvecs(log,0,"box_rel",state->box_rel,DIM);
	  if (list&LIST_V)
		  pr_rvecs(log,0,"boxv",state->boxv,DIM);
	  if (list&LIST_X)
		  pr_rvecs(log,0,int_title("x",0,buf,255),state->x,state->natoms);
	  if (list&LIST_V)
		  pr_rvecs(log,0,int_title("v",0,buf,255),state->v,state->natoms);
	  if (list&LIST_TOP)
		  pr_mtop(log,0,int_title("topology",cr->nodeid,buf,255),mtop,TRUE);
	  fflush(log);
  }
}

static void pr_molblock(FILE *fp,int indent,const char *title,
                        gmx_molblock_t *molb,int n,
                        gmx_moltype_t *molt,
                        gmx_bool bShowNumbers)
{
    indent = pr_title_n(fp,indent,title,n);
    (void) pr_indent(fp,indent);
    (void) fprintf(fp,"%-20s = %d \"%s\"\n",
                   "moltype",molb->type,*(molt[molb->type].name));
    pr_int(fp,indent,"#molecules",molb->nmol);
    pr_int(fp,indent,"#atoms_mol",molb->natoms_mol);
    pr_int(fp,indent,"#posres_xA",molb->nposres_xA);
    if (molb->nposres_xA > 0) {
        pr_rvecs(fp,indent,"posres_xA",molb->posres_xA,molb->nposres_xA);
    }
    pr_int(fp,indent,"#posres_xB",molb->nposres_xB);
    if (molb->nposres_xB > 0) {
        pr_rvecs(fp,indent,"posres_xB",molb->posres_xB,molb->nposres_xB);
    }
}

static void check_pbc(FILE *fp,rvec x[],int shell)
{
  int m,now;
  
  now = shell-4;
  for(m=0; (m<DIM); m++)
    if (fabs(x[shell][m]-x[now][m]) > 0.3) {
      pr_rvecs(fp,0,"SHELL-X",x+now,5);
      break;
    }
}

void list_xtc(char *fn, bool bXVG)
{
  int    xd,indent;
  char   buf[256];
  rvec   *x;
  matrix box;
  int    nframe,natoms,step;
  real   prec,time;
  bool   bOK;
  
  xd = open_xtc(fn,"r");
  read_first_xtc(xd,&natoms,&step,&time,box,&x,&prec,&bOK);
		
  nframe=0;
  do {
    if (bXVG) {
      int i,d;
      
      fprintf(stdout,"%g",time);
      for(i=0; i<natoms; i++)
	for(d=0; d<DIM; d++)
	  fprintf(stdout," %g",x[i][d]);
      fprintf(stdout,"\n");
    } else {
      sprintf(buf,"%s frame %d",fn,nframe);
      indent=0;
      indent=pr_title(stdout,indent,buf);
      pr_indent(stdout,indent);
      fprintf(stdout,"natoms=%10d  step=%10d  time=%12.7e  prec=%10g\n",
	    natoms,step,time,prec);
      pr_rvecs(stdout,indent,"box",box,DIM);
      pr_rvecs(stdout,indent,"x",x,natoms);
    }
    nframe++;
  } while (read_next_xtc(xd,natoms,&step,&time,box,x,&prec,&bOK));
  if (!bOK)
    fprintf(stderr,"\nWARNING: Incomplete frame at time %g\n",time);
  close_xtc(xd);
}

real calc_gyro(rvec x[], int gnx, atom_id index[], t_atom atom[], real tm,
               rvec gvec, rvec d, gmx_bool bQ, gmx_bool bRot, gmx_bool bMOI, matrix trans)
{
    int    i, ii, m;
    real   gyro, dx2, m0, Itot;
    rvec   comp;

    if (bRot)
    {
        principal_comp(gnx, index, atom, x, trans, d);
        Itot = norm(d);
        if (bMOI)
        {
            return Itot;
        }
        for (m = 0; (m < DIM); m++)
        {
            d[m] = std::sqrt(d[m]/tm);
        }
#ifdef DEBUG
        pr_rvecs(stderr, 0, "trans", trans, DIM);
#endif
        /* rotate_atoms(gnx,index,x,trans); */
    }
    clear_rvec(comp);
    for (i = 0; (i < gnx); i++)
    {
        ii = index[i];
        if (bQ)
        {
            m0 = std::abs(atom[ii].q);
        }
        else
        {
            m0 = atom[ii].m;
        }
        for (m = 0; (m < DIM); m++)
        {
            dx2      = x[ii][m]*x[ii][m];
            comp[m] += dx2*m0;
        }
    }
    gyro = comp[XX]+comp[YY]+comp[ZZ];

    for (m = 0; (m < DIM); m++)
    {
        gvec[m] = std::sqrt((gyro-comp[m])/tm);
    }

    return std::sqrt(gyro/tm);
}

void calc_virial(FILE *log,int start,int homenr,rvec x[],rvec f[],
		 tensor vir_part,tensor pme_vir,
		 t_graph *graph,matrix box,
		 t_nrnb *nrnb,t_forcerec *fr,bool bTweak)
{
  int i,j;
  tensor virtest;
  
  /* Now it is time for the short range virial. At this timepoint vir_part
   * already contains the virial from surrounding boxes.
   * Calculate partial virial, for local atoms only, based on short range. 
   * Total virial is computed in global_stat, called from do_md 
   */
  f_calc_vir(log,start,start+homenr,x,f,vir_part,graph,box);
  inc_nrnb(nrnb,eNR_VIRIAL,homenr);

  /* Add up the long range forces if necessary */
  /* if (!bTweak) {
    sum_lrforces(f,fr,start,homenr);
    }*/
  
  /* Add up virial if necessary */  
  if (EEL_LR(fr->eeltype) && (fr->eeltype != eelPPPM)) {
    if (debug && bTweak) {
      clear_mat(virtest);
      f_calc_vir(log,start,start+homenr,x,fr->f_pme,virtest,graph,box);
      pr_rvecs(debug,0,"virtest",virtest,DIM);
      pr_rvecs(debug,0,"pme_vir",pme_vir,DIM);
    }    
    /* PPPM virial sucks */
    if (!bTweak)
      for(i=0; (i<DIM); i++) 
	for(j=0; (j<DIM); j++) 
	  vir_part[i][j]+=pme_vir[i][j];
  }
  if (debug)
    pr_rvecs(debug,0,"vir_part",vir_part,DIM);
}