C++ clear_mat函数代码示例

void init_state(t_state *state,int natoms,int ngtc)
{
  int i;

  state->natoms = natoms;
  state->nrng   = 0;
  state->flags  = 0;
  state->lambda = 0;
  clear_mat(state->box);
  clear_mat(state->box_rel);
  clear_mat(state->boxv);
  clear_mat(state->pres_prev);
  init_gtc_state(state,ngtc);
  state->nalloc = state->natoms;
  if (state->nalloc > 0) {
    snew(state->x,state->nalloc);
    snew(state->v,state->nalloc);
  } else {
    state->x = NULL;
    state->v = NULL;
  }
  state->sd_X = NULL;
  state->cg_p = NULL;

  init_ekinstate(&state->ekinstate);

  init_energyhistory(&state->enerhist);

  state->ddp_count = 0;
  state->ddp_count_cg_gl = 0;
  state->cg_gl = NULL;
  state->cg_gl_nalloc = 0;
}

static void init_grptcstat(int ngtc, t_grp_tcstat tcstat[])
{
    int i, j;

    for (i = 0; (i < ngtc); i++)
    {
        tcstat[i].T = 0;
        clear_mat(tcstat[i].ekinh);
        clear_mat(tcstat[i].ekinh_old);
        clear_mat(tcstat[i].ekinf);
    }
}

static void calc_x_times_f(int nxf, const rvec x[], const rvec f[],
                           gmx_bool bScrewPBC, const matrix box,
                           matrix x_times_f)
{
    clear_mat(x_times_f);

    for (int i = 0; i < nxf; i++)
    {
        for (int d = 0; d < DIM; d++)
        {
            for (int n = 0; n < DIM; n++)
            {
                x_times_f[d][n] += x[i][d]*f[i][n];
            }
        }

        if (bScrewPBC)
        {
            int isx = IS2X(i);
            /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
            if (isx == 1 || isx == -1)
            {
                for (int d = 0; d < DIM; d++)
                {
                    for (int n = 0; n < DIM; n++)
                    {
                        x_times_f[d][n] += box[d][d]*f[i][n];
                    }
                }
            }
        }
    }
}

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 int xyz_first_x(t_trxstatus *status, FILE *fp, const output_env_t oenv,
                       real *t, rvec **x, matrix box)
/* Reads fp, mallocs x, and returns x and box
 * Returns natoms when successful, FALSE otherwise
 */
{
    int    m;

    initcount(status);

    clear_mat(box);
    choose_file_format(fp);

    for (m = 0; (m < DIM); m++)
    {
        box[m][m] = status->BOX[m];
    }

    snew(*x, status->NATOMS);
    *t = status->DT;
    if (!xyz_next_x(status, fp, oenv, t, status->NATOMS, *x, box))
    {
        return 0;
    }
    *t = 0.0;

    return status->NATOMS;
}

static void rotate_x(int natom,rvec xin[],real angle,rvec xout[],
		     gmx_bool bZ,gmx_bool bUpsideDown,real dz)
{
  int i;
  matrix mat;
  
  angle *= DEG2RAD;
  clear_mat(mat);
  if (bZ) {
    mat[XX][XX] = cos(angle);
    mat[XX][YY] = sin(angle);
    mat[YY][XX] = -sin(angle);
    mat[YY][YY] = cos(angle);
    mat[ZZ][ZZ] = 1;
  }
  else {
    mat[XX][XX] = 1;
    mat[YY][YY] = cos(angle);
    mat[YY][ZZ] = sin(angle);
    mat[ZZ][YY] = -sin(angle);
    mat[ZZ][ZZ] = cos(angle);
  }
    
  for(i=0; (i<natom); i++) {
    mvmul(mat,xin[i],xout[i]);
    if (bUpsideDown)
      xout[i][ZZ] *= -1;
    xout[i][ZZ] += dz;
  }
}

void clear_trxframe(t_trxframe *fr,gmx_bool bFirst)
{
  fr->not_ok  = 0;
  fr->bTitle  = FALSE;
  fr->bStep   = FALSE;
  fr->bTime   = FALSE;
  fr->bLambda = FALSE;
  fr->bAtoms  = FALSE;
  fr->bPrec   = FALSE;
  fr->bX      = FALSE;
  fr->bV      = FALSE;
  fr->bF      = FALSE;
  fr->bBox    = FALSE;
  if (bFirst) {
    fr->flags  = 0;
    fr->bDouble= FALSE;
    fr->natoms = -1;
    fr->t0     = 0;
    fr->tpf    = 0;
    fr->tppf   = 0;
    fr->title  = NULL;
    fr->step   = 0;
    fr->time   = 0;
    fr->lambda = 0;
    fr->atoms  = NULL;
    fr->prec   = 0;
    fr->x      = NULL;
    fr->v      = NULL;
    fr->f      = NULL;
    clear_mat(fr->box);
    fr->bPBC   = FALSE;
    fr->ePBC   = -1;
  }
}

PyObject *wrap_calc_fit_R(PyObject *self,PyObject *args)
{

  PyObject *cs1, *cs2, *mass;
  if(!PyArg_ParseTuple(args,"OOO",&cs1, &cs2, &mass))
    return NULL;
  int natoms1 = PySequence_Length(cs1);
  int natoms2 = PySequence_Length(cs2);
  if( natoms1 != natoms2 ) {
    Error("Cannot fit coordinate sets with different lengths");
  }
  rvec x1[natoms1];
  rvec x2[natoms1];
  real m[natoms1];
  PyObject2rvec( cs1, x1, natoms1);
  PyObject2rvec( cs2, x2, natoms2);
  PyObject2real_array(mass, m, natoms1);
  center(x1, natoms1);
  center(x2, natoms1);
  matrix R;
  clear_mat(R);
  calc_fit_R(natoms1,m,x1,x2,R);
  PyObject *ret = matrix2PyObject(R);
  return ret;
}

// -----------------------------------------------------------------
// Take log of hermitian part of decomposition to define scalar field
void matrix_log(matrix *in, matrix *out) {
  char V = 'V';     // Ask LAPACK for both eigenvalues and eigenvectors
  char U = 'U';     // Have LAPACK store upper triangle of in
  int row, col, Npt = NCOL, stat = 0, Nwork = 2 * NCOL;
  matrix evecs, tmat;

  // Convert in to column-major double array used by LAPACK
  for (row = 0; row < NCOL; row++) {
    for (col = 0; col < NCOL; col++) {
      store[2 * (col * NCOL + row)] = in->e[row][col].real;
      store[2 * (col * NCOL + row) + 1] = in->e[row][col].imag;
    }
  }

  // Compute eigenvalues and eigenvectors of in
  zheev_(&V, &U, &Npt, store, &Npt, eigs, work, &Nwork, Rwork, &stat);

  if (stat != 0)
    printf("WARNING: zheev returned error message %d\n", stat);

  // Move the results back into matrix structures
  // Use evecs to hold the eigenvectors for projection
  clear_mat(out);
  for (row = 0; row < NCOL; row++) {
    for (col = 0; col < NCOL; col++) {
      evecs.e[row][col].real = store[2 * (col * NCOL + row)];
      evecs.e[row][col].imag = store[2 * (col * NCOL + row) + 1];
    }
    out->e[row][row].real = log(eigs[row]);
  }
  // Inverse of eigenvector matrix is simply adjoint
  mult_na(out, &evecs, &tmat);
  mult_nn(&evecs, &tmat, out);
}

void clear_trxframe(t_trxframe *fr, gmx_bool bFirst)
{
    fr->not_ok    = 0;
    fr->bStep     = FALSE;
    fr->bTime     = FALSE;
    fr->bLambda   = FALSE;
    fr->bFepState = FALSE;
    fr->bAtoms    = FALSE;
    fr->bPrec     = FALSE;
    fr->bX        = FALSE;
    fr->bV        = FALSE;
    fr->bF        = FALSE;
    fr->bBox      = FALSE;
    if (bFirst)
    {
        fr->bDouble   = FALSE;
        fr->natoms    = -1;
        fr->step      = 0;
        fr->time      = 0;
        fr->lambda    = 0;
        fr->fep_state = 0;
        fr->atoms     = nullptr;
        fr->prec      = 0;
        fr->x         = nullptr;
        fr->v         = nullptr;
        fr->f         = nullptr;
        clear_mat(fr->box);
        fr->bPBC   = FALSE;
        fr->ePBC   = -1;
    }
}

/* Center of mass code for groups */
void calc_vcm_grp(int start, int homenr, t_mdatoms *md,
                  rvec x[], rvec v[], t_vcm *vcm)
{
    int    i, g, m;
    real   m0;
    rvec   j0;

    if (vcm->mode != ecmNO)
    {
        /* Also clear a possible rest group */
        for (g = 0; (g < vcm->nr+1); g++)
        {
            /* Reset linear momentum */
            vcm->group_mass[g] = 0;
            clear_rvec(vcm->group_p[g]);

            if (vcm->mode == ecmANGULAR)
            {
                /* Reset anular momentum */
                clear_rvec(vcm->group_j[g]);
                clear_rvec(vcm->group_x[g]);
                clear_rvec(vcm->group_w[g]);
                clear_mat(vcm->group_i[g]);
            }
        }

        g = 0;
        for (i = start; (i < start+homenr); i++)
        {
            m0 = md->massT[i];
            if (md->cVCM)
            {
                g = md->cVCM[i];
            }

            /* Calculate linear momentum */
            vcm->group_mass[g]  += m0;
            for (m = 0; (m < DIM); m++)
            {
                vcm->group_p[g][m] += m0*v[i][m];
            }

            if (vcm->mode == ecmANGULAR)
            {
                /* Calculate angular momentum */
                cprod(x[i], v[i], j0);

                for (m = 0; (m < DIM); m++)
                {
                    vcm->group_j[g][m] += m0*j0[m];
                    vcm->group_x[g][m] += m0*x[i][m];
                }
                /* Update inertia tensor */
                update_tensor(x[i], m0, vcm->group_i[g]);
            }
        }
    }
}

void get_enx_state(char *fn, real t, gmx_groups_t *groups, t_inputrec *ir,
                   t_state *state)
{
  /* Should match the names in mdebin.c */
  static const char *boxvel_nm[] = {
  "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
  "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
  };
  
  static const char *pcouplmu_nm[] = {
    "Pcoupl-Mu-XX", "Pcoupl-Mu-YY", "Pcoupl-Mu-ZZ",
    "Pcoupl-Mu-YX", "Pcoupl-Mu-ZX", "Pcoupl-Mu-ZY"
  };
  int ind0[] = { XX,YY,ZZ,YY,ZZ,ZZ };
  int ind1[] = { XX,YY,ZZ,XX,XX,YY };

  int in,nre,nfr,i,ni,npcoupl;
  char       buf[STRLEN];
  gmx_enxnm_t *enm;
  t_enxframe *fr;

  in = open_enx(fn,"r");
  do_enxnms(in,&nre,&enm);
  snew(fr,1);
  nfr = 0;
  while ((nfr==0 || fr->t != t) && do_enx(in,fr)) {
    nfr++;
  }
  close_enx(in);
  fprintf(stderr,"\n");

  if (nfr == 0 || fr->t != t)
    gmx_fatal(FARGS,"Could not find frame with time %f in '%s'",t,fn);
  
  npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
  if (ir->epc == epcPARRINELLORAHMAN) {
    clear_mat(state->boxv);
    for(i=0; i<npcoupl; i++) {
      state->boxv[ind0[i]][ind1[i]] =
	find_energy(boxvel_nm[i],nre,enm,fr);
    }
    fprintf(stderr,"\nREAD %d BOX VELOCITIES FROM %s\n\n",npcoupl,fn);
  }

  if (ir->etc == etcNOSEHOOVER) {
    for(i=0; i<state->ngtc; i++) {
      ni = groups->grps[egcTC].nm_ind[i];
      sprintf(buf,"Xi-%s",*(groups->grpname[ni]));
      state->nosehoover_xi[i] = find_energy(buf,nre,enm,fr);
    }
    fprintf(stderr,"\nREAD %d NOSE-HOOVER Xi's FROM %s\n\n",state->ngtc,fn);
  }
  
  free_enxnms(nre,enm);
  free_enxframe(fr);
  sfree(fr);
}

real calc_cm(FILE *log,int natoms,real mass[],rvec x[],rvec v[],
	     rvec xcm,rvec vcm,rvec acm,matrix L)
{
  rvec dx,a0;
  real tm,m0;
  int  i,m;

  clear_rvec(xcm);
  clear_rvec(vcm);
  clear_rvec(acm);
  tm=0.0;
  for(i=0; (i<natoms); i++) {
    m0=mass[i];
    tm+=m0;
    cprod(x[i],v[i],a0);
    for(m=0; (m<DIM); m++) {
      xcm[m]+=m0*x[i][m]; /* c.o.m. position */
      vcm[m]+=m0*v[i][m]; /* c.o.m. velocity */
      acm[m]+=m0*a0[m];   /* rotational velocity around c.o.m. */
    }
  }
  cprod(xcm,vcm,a0);
  for(m=0; (m<DIM); m++) {
    xcm[m]/=tm;
    vcm[m]/=tm;
    acm[m]-=a0[m]/tm;
  }

#define PVEC(str,v) fprintf(log,\
			    "%s[X]: %10.5e  %s[Y]: %10.5e  %s[Z]: %10.5e\n", \
			   str,v[0],str,v[1],str,v[2])
#ifdef DEBUG
  PVEC("xcm",xcm);
  PVEC("acm",acm);
  PVEC("vcm",vcm);
#endif
  
  clear_mat(L);
  for(i=0; (i<natoms); i++) {
    m0=mass[i];
    for(m=0; (m<DIM); m++)
      dx[m]=x[i][m]-xcm[m];
    L[XX][XX]+=dx[XX]*dx[XX]*m0;
    L[XX][YY]+=dx[XX]*dx[YY]*m0;
    L[XX][ZZ]+=dx[XX]*dx[ZZ]*m0;
    L[YY][YY]+=dx[YY]*dx[YY]*m0;
    L[YY][ZZ]+=dx[YY]*dx[ZZ]*m0;
    L[ZZ][ZZ]+=dx[ZZ]*dx[ZZ]*m0;
  }
#ifdef DEBUG
  PVEC("L-x",L[XX]);
  PVEC("L-y",L[YY]);
  PVEC("L-z",L[ZZ]);
#endif

  return tm;
}

/* the non-mass-weighted mean-squared displacement calcuation */
static real calc1_norm(t_corr *curr, int nx, atom_id index[], int nx0, rvec xc[],
                       rvec dcom, gmx_bool bTen, matrix mat)
{
    int  i, ix, m, m2;
    real g, r, r2;
    rvec rv;

    g = 0.0;
    clear_mat(mat);

    for (i = 0; (i < nx); i++)
    {
        ix = index[i];
        r2 = 0.0;
        switch (curr->type)
        {
            case NORMAL:
                for (m = 0; (m < DIM); m++)
                {
                    rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
                    r2   += rv[m]*rv[m];
                    if (bTen)
                    {
                        for (m2 = 0; m2 <= m; m2++)
                        {
                            mat[m][m2] += rv[m]*rv[m2];
                        }
                    }
                }
                break;
            case X:
            case Y:
            case Z:
                r = xc[ix][curr->type-X] - curr->x0[nx0][ix][curr->type-X] -
                    dcom[curr->type-X];
                r2 += r*r;
                break;
            case LATERAL:
                for (m = 0; (m < DIM); m++)
                {
                    if (m != curr->axis)
                    {
                        r   = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
                        r2 += r*r;
                    }
                }
                break;
            default:
                gmx_fatal(FARGS, "Error: did not expect option value %d", curr->type);
        }
        g += r2;
    }
    g /= nx;
    msmul(mat, 1.0/nx, mat);

    return g;
}

void set_box_rel(t_inputrec *ir,t_state *state)
{
  /* Make sure the box obeys the restrictions before we fix the ratios */
  correct_box(NULL,0,state->box,NULL);

  clear_mat(state->box_rel);

  if (PRESERVE_SHAPE(*ir))
    do_box_rel(ir,state->box_rel,state->box,TRUE);
}