static void calc_axes(rvec x[],t_atom atom[],int gnx[],atom_id *index[],
		      gmx_bool bRot,t_bundle *bun)
{
  int  end,i,div,d;
  real *mtot,m;
  rvec axis[MAX_ENDS],cent;
  
  snew(mtot,bun->n);

  for(end=0; end<bun->nend; end++) {
    for(i=0; i<bun->n; i++) {
      clear_rvec(bun->end[end][i]);
      mtot[i] = 0;
    }
    div = gnx[end]/bun->n;
    for(i=0; i<gnx[end]; i++) {
      m = atom[index[end][i]].m;
      for(d=0; d<DIM; d++)
	bun->end[end][i/div][d] += m*x[index[end][i]][d];
      mtot[i/div] += m;
    }
    clear_rvec(axis[end]);
    for(i=0; i<bun->n; i++) {
      svmul(1.0/mtot[i],bun->end[end][i],bun->end[end][i]);
      rvec_inc(axis[end],bun->end[end][i]);
    }
    svmul(1.0/bun->n,axis[end],axis[end]);
  }
  sfree(mtot);

  rvec_add(axis[0],axis[1],cent);
  svmul(0.5,cent,cent);
  /* center the bundle on the origin */
  for(end=0; end<bun->nend; end++) {
    rvec_dec(axis[end],cent);
    for(i=0; i<bun->n; i++)
      rvec_dec(bun->end[end][i],cent);
  }
  if (bRot) {
    /* rotate the axis parallel to the z-axis */
    rotate_ends(bun,axis[0],YY,ZZ);
    rotate_ends(bun,axis[0],XX,ZZ);
  }
  for(i=0; i<bun->n; i++) {
    rvec_add(bun->end[0][i],bun->end[1][i],bun->mid[i]);
    svmul(0.5,bun->mid[i],bun->mid[i]);
    rvec_sub(bun->end[0][i],bun->end[1][i],bun->dir[i]);
    bun->len[i] = norm(bun->dir[i]);
    unitv(bun->dir[i],bun->dir[i]);
  }
}

void calc_triclinic_images(matrix box,rvec img[])
{
    int i;

    /* Calculate 3 adjacent images in the xy-plane */
    copy_rvec(box[0],img[0]);
    copy_rvec(box[1],img[1]);
    if (img[1][XX] < 0)
        svmul(-1,img[1],img[1]);
    rvec_sub(img[1],img[0],img[2]);

    /* Get the next 3 in the xy-plane as mirror images */
    for(i=0; i<3; i++)
        svmul(-1,img[i],img[3+i]);

    /* Calculate the first 4 out of xy-plane images */
    copy_rvec(box[2],img[6]);
    if (img[6][XX] < 0)
        svmul(-1,img[6],img[6]);
    for(i=0; i<3; i++)
        rvec_add(img[6],img[i+1],img[7+i]);

    /* Mirror the last 4 from the previous in opposite rotation */
    for(i=0; i<4; i++)
        svmul(-1,img[6 + (2+i) % 4],img[10+i]);
}

PyObject *wrap_fit(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);
  rvec cent;
  center_and_get_vec(x1, natoms1, cent);     // center x1 and get vector for back translation
  center(x2, natoms1);                                // center x2
  do_fit(natoms1, m, x1, x2);               

  int i;
  for(i=0;i<natoms1;i++)                    // translate back
    rvec_add( x2[i], cent, x2[i]);

  PyObject *ret = rvec2PyObject(x2, natoms1);
  return ret;
}

static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
                          real *rmin, real *rmax, int *min_ind)
{
#define NSHIFT 26
    int  sx, sy, sz, i, j, s;
    real sqr_box, r2min, r2max, r2;
    rvec shift[NSHIFT], d0, d;

    sqr_box = sqr(min(norm(box[XX]), min(norm(box[YY]), norm(box[ZZ]))));

    s = 0;
    for (sz = -1; sz <= 1; sz++)
    {
        for (sy = -1; sy <= 1; sy++)
        {
            for (sx = -1; sx <= 1; sx++)
            {
                if (sx != 0 || sy != 0 || sz != 0)
                {
                    for (i = 0; i < DIM; i++)
                    {
                        shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
                    }
                    s++;
                }
            }
        }
    }

    r2min = sqr_box;
    r2max = 0;

    for (i = 0; i < n; i++)
    {
        for (j = i+1; j < n; j++)
        {
            rvec_sub(x[index[i]], x[index[j]], d0);
            r2 = norm2(d0);
            if (r2 > r2max)
            {
                r2max = r2;
            }
            for (s = 0; s < NSHIFT; s++)
            {
                rvec_add(d0, shift[s], d);
                r2 = norm2(d);
                if (r2 < r2min)
                {
                    r2min      = r2;
                    min_ind[0] = i;
                    min_ind[1] = j;
                }
            }
        }
    }

    *rmin = sqrt(r2min);
    *rmax = sqrt(r2max);
}

void visualize_box(FILE *out,int a0,int r0,matrix box,rvec gridsize)
{
  int     *edge;
  rvec    *vert,shift;
  int     nx,ny,nz,nbox,nat;
  int     i,j,x,y,z;
  int     rectedge[24] = { 0,1, 1,3, 3,2, 0,2, 0,4, 1,5, 3,7, 2,6, 4,5, 5,7, 7,6, 6,4 };

  a0++;
  r0++;
  
  nx = (int)(gridsize[XX]+0.5);
  ny = (int)(gridsize[YY]+0.5);
  nz = (int)(gridsize[ZZ]+0.5);
  nbox = nx*ny*nz;
  if (TRICLINIC(box)) {
    nat = nbox*NCUCVERT;
    snew(vert,nat);
    calc_compact_unitcell_vertices(ecenterDEF,box,vert);
    j = 0;
    for(z=0; z<nz; z++)
      for(y=0; y<ny; y++)
	for(x=0; x<nx; x++) {
	  for(i=0; i<DIM; i++)
	    shift[i] = x*box[0][i]+y*box[1][i]+z*box[2][i];
	  for(i=0; i<NCUCVERT; i++) {
	    rvec_add(vert[i],shift,vert[j]);
	    j++;
	  }
	}
    
    for(i=0; i<nat; i++) {
      fprintf(out,pdbformat,"ATOM",a0+i,"C","BOX",'K'+i/NCUCVERT,r0+i,
	      10*vert[i][XX],10*vert[i][YY],10*vert[i][ZZ]);
      fprintf(out,"\n");
    }
    
    edge = compact_unitcell_edges();
    for(j=0; j<nbox; j++)
      for(i=0; i<NCUCEDGE; i++)
	fprintf(out,"CONECT%5d%5d\n",
		a0 + j*NCUCVERT + edge[2*i],
		a0 + j*NCUCVERT + edge[2*i+1]);
    
    sfree(vert);
  } else {
    i=0;
    for(z=0; z<=1; z++)
      for(y=0; y<=1; y++)
	for(x=0; x<=1; x++) {
	  fprintf(out,pdbformat,"ATOM",a0+i,"C","BOX",'K'+i/8,r0+i,
		  x*10*box[XX][XX],y*10*box[YY][YY],z*10*box[ZZ][ZZ]);
	  fprintf(out,"\n");
	  i++;
	}
    for(i=0; i<24; i+=2)
      fprintf(out,"CONECT%5d%5d\n",a0+rectedge[i],a0+rectedge[i+1]);
  }
}

static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,
                         atom_id index[], rvec xref, int ePBC)
{
    const real tol = 1e-4;
    gmx_bool   bChanged;
    int        m, j, ai, iter;
    real       mass, mtot;
    rvec       dx, xtest;

    /* First simple calculation */
    clear_rvec(xref);
    mtot = 0;
    for (m = 0; (m < nrefat); m++)
    {
        ai   = index[m];
        mass = top->atoms.atom[ai].m;
        for (j = 0; (j < DIM); j++)
        {
            xref[j] += mass*x[ai][j];
        }
        mtot += mass;
    }
    svmul(1/mtot, xref, xref);
    /* Now check if any atom is more than half the box from the COM */
    if (ePBC != epbcNONE)
    {
        iter = 0;
        do
        {
            bChanged = FALSE;
            for (m = 0; (m < nrefat); m++)
            {
                ai   = index[m];
                mass = top->atoms.atom[ai].m/mtot;
                pbc_dx(pbc, x[ai], xref, dx);
                rvec_add(xref, dx, xtest);
                for (j = 0; (j < DIM); j++)
                {
                    if (std::abs(xtest[j]-x[ai][j]) > tol)
                    {
                        /* Here we have used the wrong image for contributing to the COM */
                        xref[j] += mass*(xtest[j]-x[ai][j]);
                        x[ai][j] = xtest[j];
                        bChanged = TRUE;
                    }
                }
            }
            if (bChanged)
            {
                printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
            }
            iter++;
        }
        while (bChanged);
    }
}

PyObject *apply_rotation( PyObject *self, PyObject *args)
{
  PyObject *Rotation;
  PyObject *py_v;
  real phi;
  if(!PyArg_ParseTuple(args,"OOd",&Rotation,&py_v, &phi))
    return NULL;

  PyObject *rm1 = PyObject_GetAttrString(Rotation,"m1");
  PyObject *rm2 = PyObject_GetAttrString(Rotation,"m2");

  matrix m1, m2;
  PyObject2matrix(rm1, m1);
  PyObject2matrix(rm2, m2);
  rvec v, v2;

  PyObject *py_v2 = PyObject_GetAttrString(Rotation,"v2");

  Pyvec2rvec(py_v, v);
  Pyvec2rvec(py_v2, v2);

  rvec vec;
  rvec_sub(v, v2, vec );

  rvec b, d, a, c, e;
  mvmul(m1, vec, b);
  mvmul(m2, vec, d);
  real cc = cos(phi);
  svmul( cc, vec, a);
  svmul( -cc, b, c);
  svmul( sin(phi), d, e);

  clear_rvec( vec );
  rvec_add( a, b, vec);
  rvec_add( vec, c, vec );
  rvec_add( vec, e, vec );
  
  clear_rvec(v);
  rvec_add( v2, vec, v);
  return Py_BuildValue("[ddd]", v[XX], v[YY], v[ZZ] );
 }

void write_conf_to_grofile(const System& system,
                           const std::string& path,
                           const real sigma,
                           const OutputMode mode)
{
    constexpr char ATOM_NAME[2] = "C";
    constexpr char RESIDUE_NAME[4] = "SOL";

    auto open_mode = std::ios::out;

    if (mode == OutputMode::Append)
    {
        open_mode = open_mode | std::ios::app;
    }

    std::ofstream out { path, open_mode };

    out << system.title << '\n'
        << system.num_atoms() << '\n';

    out.setf(std::ios::fixed);
    out.precision(3);

    uint64_t n = 1;

    for (const auto &list : system.cell_lists)
    {
        for (unsigned i = 0; i < list.num_atoms(); ++i)
        {
            const auto x0 = RVec {
                list.xs.at(i * NDIM + XX),
                list.xs.at(i * NDIM + YY),
                list.xs.at(i * NDIM + ZZ)
            };
            const auto xabs = rvec_add(x0, list.origin);

            out << std::setw(5) << std::right << n
                << std::setw(5) << std::left << RESIDUE_NAME
                << std::setw(5) << std::right << ATOM_NAME
                << std::setw(5) << n
                << std::setw(8) << xabs[XX] * sigma
                << std::setw(8) << xabs[YY] * sigma
                << std::setw(8) << xabs[ZZ] * sigma
                << '\n';

            ++n;
        }
    }

    out << std::setw(9) << std::right << system.box_size[0] * sigma << ' '
        << std::setw(9) << system.box_size[1] * sigma << ' '
        << std::setw(9) << system.box_size[2] * sigma << '\n';
}

//! Helper method to calculate a vector from two or three positions.
static void
calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
{
    switch (natoms)
    {
        case 2:
            if (pbc)
            {
                pbc_dx(pbc, x[1], x[0], xout);
            }
            else
            {
                rvec_sub(x[1], x[0], xout);
            }
            svmul(0.5, xout, cout);
            rvec_add(x[0], cout, cout);
            break;
        case 3:
        {
            rvec v1, v2;
            if (pbc)
            {
                pbc_dx(pbc, x[1], x[0], v1);
                pbc_dx(pbc, x[2], x[0], v2);
            }
            else
            {
                rvec_sub(x[1], x[0], v1);
                rvec_sub(x[2], x[0], v2);
            }
            cprod(v1, v2, xout);
            rvec_add(x[0], x[1], cout);
            rvec_add(cout, x[2], cout);
            svmul(1.0/3.0, cout, cout);
            break;
        }
        default:
            GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
    }
}

/*!
 * \param[in]  top    Topology structure with masses.
 * \param[in]  x      Position vectors of all atoms.
 * \param[in]  pbc    Periodic boundary conditions structure.
 * \param[in]  nrefat Number of atoms in the index.
 * \param[in]  index  Indices of atoms.
 * \param[out] xout   COM position for the indexed atoms.
 *
 * Works as gmx_calc_com(), but takes into account periodic boundary
 * conditions: If any atom is more than half the box from the COM,
 * it is wrapped around and a new COM is calculated. This is repeated
 * until no atoms violate the condition.
 *
 * Modified from src/tools/gmx_sorient.c in Gromacs distribution.
 */
void
gmx_calc_com_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,
                 int nrefat, const int index[], rvec xout)
{
    GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
                       "No masses available while mass weighting was requested");
    /* First simple calculation */
    clear_rvec(xout);
    real mtot = 0;
    int  molb = 0;
    for (int m = 0; m < nrefat; ++m)
    {
        const int  ai   = index[m];
        const real mass = mtopGetAtomMass(top, ai, &molb);
        for (int j = 0; j < DIM; ++j)
        {
            xout[j] += mass * x[ai][j];
        }
        mtot += mass;
    }
    svmul(1.0/mtot, xout, xout);
    /* Now check if any atom is more than half the box from the COM */
    if (pbc)
    {
        const real tol  = 1e-4;
        bool       bChanged;
        do
        {
            bChanged = false;
            molb     = 0;
            for (int m = 0; m < nrefat; ++m)
            {
                rvec       dx, xtest;
                const int  ai   = index[m];
                const real mass = mtopGetAtomMass(top, ai, &molb) / mtot;
                pbc_dx(pbc, x[ai], xout, dx);
                rvec_add(xout, dx, xtest);
                for (int j = 0; j < DIM; ++j)
                {
                    if (fabs(xtest[j] - x[ai][j]) > tol)
                    {
                        /* Here we have used the wrong image for contributing to the COM */
                        xout[j] += mass * (xtest[j] - x[ai][j]);
                        x[ai][j] = xtest[j];
                        bChanged = true;
                    }
                }
            }
        }
        while (bChanged);
    }
}

static void calc_ringh(rvec xattach, rvec xb, rvec xc, rvec xh)
{
    rvec tab, tac;
    real n;

    /* Add a proton on a ring to atom attach at distance 0.1 nm */
    rvec_sub(xattach, xb, tab);
    rvec_sub(xattach, xc, tac);
    rvec_add(tab, tac, xh);
    n = 0.1/norm(xh);
    svmul(n, xh, xh);
    rvec_inc(xh, xattach);
}

static void
calc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout)
{
    switch (natoms)
    {
        case 2:
            if (pbc)
            {
                pbc_dx(pbc, x[1], x[0], xout);
            }
            else
            {
                rvec_sub(x[1], x[0], xout);
            }
            svmul(0.5, xout, cout);
            rvec_add(x[0], cout, cout);
            break;
        case 3: {
            rvec v1, v2;
            if (pbc)
            {
                pbc_dx(pbc, x[1], x[0], v1);
                pbc_dx(pbc, x[2], x[0], v2);
            }
            else
            {
                rvec_sub(x[1], x[0], v1);
                rvec_sub(x[2], x[0], v2);
            }
            cprod(v1, v2, xout);
            rvec_add(x[0], x[1], cout);
            rvec_add(cout, x[2], cout);
            svmul(1.0/3.0, cout, cout);
            break;
        }
    }
}

void dump_pbc(FILE *fp, t_pbc *pbc)
{
    rvec sum_box;

    fprintf(fp, "ePBCDX = %d\n", pbc->ePBCDX);
    pr_rvecs(fp, 0, "box", pbc->box, DIM);
    pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
    pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
    pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
    rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
    pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
    fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
    fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
    if (pbc->ntric_vec > 0)
    {
        pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE);
        pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
    }
}

const char *
put_atoms_in_compact_unitcell(int ePBC,int ecenter,matrix box,
                              int natoms,rvec x[])
                              {
    t_pbc pbc;
    rvec box_center,dx;
    int  i;

    set_pbc(&pbc,ePBC,box);
    calc_box_center(ecenter,box,box_center);
    for(i=0; i<natoms; i++) {
        pbc_dx(&pbc,x[i],box_center,dx);
        rvec_add(box_center,dx,x[i]);
    }

    return pbc.bLimitDistance ?
        "WARNING: Could not put all atoms in the compact unitcell\n"
        : NULL;
                              }

void put_atoms_in_compact_unitcell(int ePBC, int ecenter, const matrix box,
                                   int natoms, rvec x[])
{
    t_pbc pbc;
    rvec  box_center, dx;
    int   i;

    set_pbc(&pbc, ePBC, box);

    if (pbc.ePBCDX == epbcdxUNSUPPORTED)
    {
        gmx_fatal(FARGS, "Can not put atoms in compact unitcell with unsupported PBC");
    }

    calc_box_center(ecenter, box, box_center);
    for (i = 0; i < natoms; i++)
    {
        pbc_dx(&pbc, x[i], box_center, dx);
        rvec_add(box_center, dx, x[i]);
    }
}

/*!
 * \param[in]  top    Topology structure (unused, can be NULL).
 * \param[in]  x      Position vectors of all atoms.
 * \param[in]  pbc    Periodic boundary conditions structure.
 * \param[in]  nrefat Number of atoms in the index.
 * \param[in]  index  Indices of atoms.
 * \param[out] xout   COG position for the indexed atoms.
 *
 * Works exactly as gmx_calc_com_pbc(), but calculates the center of geometry.
 */
void
gmx_calc_cog_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,
                 int nrefat, const int index[], rvec xout)
{
    const real          tol = 1e-4;
    bool                bChanged;
    int                 m, j, ai, iter;
    rvec                dx, xtest;

    /* First simple calculation */
    gmx_calc_cog(top, x, nrefat, index, xout);
    /* Now check if any atom is more than half the box from the COM */
    if (pbc)
    {
        iter = 0;
        do
        {
            bChanged = false;
            for (m = 0; m < nrefat; ++m)
            {
                ai = index[m];
                pbc_dx(pbc, x[ai], xout, dx);
                rvec_add(xout, dx, xtest);
                for (j = 0; j < DIM; ++j)
                {
                    if (fabs(xtest[j] - x[ai][j]) > tol)
                    {
                        /* Here we have used the wrong image for contributing to the COM */
                        xout[j] += (xtest[j] - x[ai][j]) / nrefat;
                        x[ai][j] = xtest[j];
                        bChanged = true;
                    }
                }
            }
            iter++;
        }
        while (bChanged);
    }
}

void random_h_coords(int natmol,int nmol,rvec x[],rvec box,
		     gmx_bool bYaw,real odist,real hdist)
{
#define cx  0.81649658
#define cy  0.47140452
#define cy2 0.94280904
#define cz  0.33333333
  
  rvec xx[24] = {
    { 0,   0,         0 }, /* O1 */
    { 0,   0,         1 }, /* H relative to Oxygen */
    { cx, cy,       -cz },
    { cx, cy,       -cz }, /* O2 */
    { 0, 0,       -1    }, /* H relative to Oxygen */
    { cx,-cy,       +cz },
    { cx, cy+cy2,     0 }, /* O3 */
    { -cx, cy,    -cz   }, /* H relative to Oxygen */
    { 0,   -cy2,    -cz },
    { 0,  2*cy+cy2, -cz }, /* O4 */
    {-cx,-cy,       +cz }, /* H relative to Oxygen */
    { 0 , cy2,      +cz },
    { 0,   0,         1 }, /* O5 */
    {-cx, cy,       +cz }, /* H relative to Oxygen */
    { 0 , -cy2,     +cz },
    { cx, cy,      1+cz }, /* O6 */
    { -cx, -cy,     -cz }, /* H relative to Oxygen */
    { 0,   cy2,     -cz },
    { cx, cy+cy2,     1 }, /* O7 */
    { 0,  0,       -1   }, /* H relative to Oxygen */
    { cx, cy,       +cz },
    { 0,  2*cy+cy2,1+cz }, /* O8 */
    { 0,   0,         1 }, /* H relative to Oxygen */
    { cx,   -cy,    -cz }
  };
  int  i,iin,iout,j,m;
  rvec tmp,t2,dip;
  
  clear_rvec(dip);
  for(i=0; (i<nmol); i++) {
    iin = natmol*i;
    iout = iin;
    svmul(odist,x[iin],x[iout]);
    svmul(-0.82,x[iout],t2);
    rvec_inc(dip,t2);
    for(j=1; (j<=2); j++) {
      svmul(hdist,xx[3*(i % 8)+j],tmp);
      rvec_add(x[iout],tmp,x[iout+j]);
      svmul(0.41,x[iout+j],t2);
      rvec_inc(dip,t2);
    }
  }
  
  box[XX] = 2*cx;
  box[YY] = 2*(cy2+cy);
  box[ZZ] = 2*(1+cz);
  for(i=0; (i<DIM); i++)
    box[i] *= odist;
    
  printf("Unitcell:  %10.5f  %10.5f  %10.5f\n",box[XX],box[YY],box[ZZ]);
  printf("Dipole:    %10.5f  %10.5f  %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
}

gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (
        gmx_sans_t  *gsans,
        rvec        *x,
        matrix       box,
        atom_id     *index,
        int          isize,
        double       binwidth,
        gmx_bool     bMC,
        gmx_bool     bNORM,
        real         mcover,
        unsigned int seed)
{
    gmx_radial_distribution_histogram_t    *pr = NULL;
    rvec              dist;
    double            rmax;
    int               i, j;
#ifdef GMX_OPENMP
    double          **tgr;
    int               tid;
    int               nthreads;
    gmx_rng_t        *trng = NULL;
#endif
    gmx_large_int_t   mc  = 0, max;
    gmx_rng_t         rng = NULL;

    /* allocate memory for pr */
    snew(pr, 1);
    /* set some fields */
    pr->binwidth = binwidth;

    /*
     * create max dist rvec
     * dist = box[xx] + box[yy] + box[zz]
     */
    rvec_add(box[XX], box[YY], dist);
    rvec_add(box[ZZ], dist, dist);

    rmax = norm(dist);

    pr->grn = (int)floor(rmax/pr->binwidth)+1;
    rmax    = pr->grn*pr->binwidth;

    snew(pr->gr, pr->grn);

    if (bMC)
    {
        /* Special case for setting automaticaly number of mc iterations to 1% of total number of direct iterations */
        if (mcover == -1)
        {
            max = (gmx_large_int_t)floor(0.5*0.01*isize*(isize-1));
        }
        else
        {
            max = (gmx_large_int_t)floor(0.5*mcover*isize*(isize-1));
        }
        rng = gmx_rng_init(seed);
#ifdef GMX_OPENMP
        nthreads = gmx_omp_get_max_threads();
        snew(tgr, nthreads);
        snew(trng, nthreads);
        for (i = 0; i < nthreads; i++)
        {
            snew(tgr[i], pr->grn);
            trng[i] = gmx_rng_init(gmx_rng_uniform_uint32(rng));
        }
        #pragma omp parallel shared(tgr,trng,mc) private(tid,i,j)
        {
            tid = gmx_omp_get_thread_num();
            /* now starting parallel threads */
            #pragma omp for
            for (mc = 0; mc < max; mc++)
            {
                i = (int)floor(gmx_rng_uniform_real(trng[tid])*isize);
                j = (int)floor(gmx_rng_uniform_real(trng[tid])*isize);
                if (i != j)
                {
                    tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
                }
            }
        }
        /* collecting data from threads */
        for (i = 0; i < pr->grn; i++)
        {
            for (j = 0; j < nthreads; j++)
            {
                pr->gr[i] += tgr[j][i];
            }
        }
        /* freeing memory for tgr and destroying trng */
        for (i = 0; i < nthreads; i++)
        {
            sfree(tgr[i]);
            gmx_rng_destroy(trng[i]);
        }
        sfree(tgr);
        sfree(trng);
#else
        for (mc = 0; mc < max; mc++)
        {
            i = (int)floor(gmx_rng_uniform_real(rng)*isize);
//.........这里部分代码省略.........

void add_conf(t_atoms *atoms, rvec **x, rvec **v, real **r, gmx_bool bSrenew,
              int ePBC, matrix box, gmx_bool bInsert,
              t_atoms *atoms_solvt, rvec *x_solvt, rvec *v_solvt, real *r_solvt,
              gmx_bool bVerbose, real rshell, int max_sol, const output_env_t oenv)
{
    t_nblist       *nlist;
    t_atoms        *atoms_all;
    real            max_vdw, *r_prot, *r_all, n2, r2, ib1, ib2;
    int             natoms_prot, natoms_solvt;
    int             i, j, jj, m, j0, j1, jjj, jnres, jnr, inr, iprot, is1, is2;
    int             prev, resnr, nresadd, d, k, ncells, maxincell;
    int             dx0, dx1, dy0, dy1, dz0, dz1;
    int             ntest, nremove, nkeep;
    rvec            dx, xi, xj, xpp, *x_all, *v_all;
    gmx_bool       *remove, *keep;
    int             bSolSol;

    natoms_prot  = atoms->nr;
    natoms_solvt = atoms_solvt->nr;
    if (natoms_solvt <= 0)
    {
        fprintf(stderr, "WARNING: Nothing to add\n");
        return;
    }

    if (ePBC == epbcSCREW)
    {
        gmx_fatal(FARGS, "Sorry, %s pbc is not yet supported", epbc_names[ePBC]);
    }

    if (bVerbose)
    {
        fprintf(stderr, "Calculating Overlap...\n");
    }

    /* Set margin around box edges to largest solvent dimension.
     * The maximum distance between atoms in a solvent molecule should
     * be calculated. At the moment a fudge factor of 3 is used.
     */
    r_prot     = *r;
    box_margin = 3*find_max_real(natoms_solvt, r_solvt);
    max_vdw    = max(3*find_max_real(natoms_prot, r_prot), box_margin);
    fprintf(stderr, "box_margin = %g\n", box_margin);

    snew(remove, natoms_solvt);

    nremove = 0;
    if (!bInsert)
    {
        for (i = 0; i < atoms_solvt->nr; i++)
        {
            if (outside_box_plus_margin(x_solvt[i], box) )
            {
                i = mark_res(i, remove, atoms_solvt->nr, atoms_solvt->atom, &nremove);
            }
        }
        fprintf(stderr, "Removed %d atoms that were outside the box\n", nremove);
    }

    /* Define grid stuff for genbox */
    /* Largest VDW radius */
    snew(r_all, natoms_prot+natoms_solvt);
    for (i = j = 0; i < natoms_prot; i++, j++)
    {
        r_all[j] = r_prot[i];
    }
    for (i = 0; i < natoms_solvt; i++, j++)
    {
        r_all[j] = r_solvt[i];
    }

    /* Combine arrays */
    combine_atoms(atoms, atoms_solvt, *x, v ? *v : NULL, x_solvt, v_solvt,
                  &atoms_all, &x_all, &v_all);

    /* Do neighboursearching step */
    do_nsgrid(stdout, bVerbose, box, x_all, atoms_all, max_vdw, oenv);

    /* check solvent with solute */
    nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
    fprintf(stderr, "nri = %d, nrj = %d\n", nlist->nri, nlist->nrj);
    for (bSolSol = 0; (bSolSol <= (bInsert ? 0 : 1)); bSolSol++)
    {
        ntest = nremove = 0;
        fprintf(stderr, "Checking %s-Solvent overlap:",
                bSolSol ? "Solvent" : "Protein");
        for (i = 0; (i < nlist->nri && nremove < natoms_solvt); i++)
        {
            inr = nlist->iinr[i];
            j0  = nlist->jindex[i];
            j1  = nlist->jindex[i+1];
            rvec_add(x_all[inr], fr->shift_vec[nlist->shift[i]], xi);

            for (j = j0; (j < j1 && nremove < natoms_solvt); j++)
            {
                jnr = nlist->jjnr[j];
                copy_rvec(x_all[jnr], xj);

                /* Check solvent-protein and solvent-solvent */
                is1 = inr-natoms_prot;
//.........这里部分代码省略.........

void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
                      matrix box,
                      rvec *x0,
                      rvec *x1, gmx_bool bX1IsCoord)
{
    gmx_specatsend_t *spas;
    rvec             *x, *vbuf, *rbuf;
    int               nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
    gmx_bool          bPBC, bScrew = FALSE;
    rvec              shift = {0, 0, 0};

    nvec = 1;
    if (x1 != NULL)
    {
        nvec++;
    }

    n = spac->at_start;
    for (d = 0; d < dd->ndim; d++)
    {
        dim = dd->dim[d];
        if (dd->nc[dim] > 2)
        {
            /* Pulse the grid forward and backward */
            vbuf = spac->vbuf;
            for (dir = 0; dir < 2; dir++)
            {
                if (dir == 0 && dd->ci[dim] == 0)
                {
                    bPBC   = TRUE;
                    bScrew = (dd->bScrewPBC && dim == XX);
                    copy_rvec(box[dim], shift);
                }
                else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
                {
                    bPBC   = TRUE;
                    bScrew = (dd->bScrewPBC && dim == XX);
                    for (i = 0; i < DIM; i++)
                    {
                        shift[i] = -box[dim][i];
                    }
                }
                else
                {
                    bPBC   = FALSE;
                    bScrew = FALSE;
                }
                spas = &spac->spas[d][dir];
                for (v = 0; v < nvec; v++)
                {
                    x = (v == 0 ? x0 : x1);
                    /* Copy the required coordinates to the send buffer */
                    if (!bPBC || (v == 1 && !bX1IsCoord))
                    {
                        /* Only copy */
                        for (i = 0; i < spas->nsend; i++)
                        {
                            copy_rvec(x[spas->a[i]], *vbuf);
                            vbuf++;
                        }
                    }
                    else if (!bScrew)
                    {
                        /* Shift coordinates */
                        for (i = 0; i < spas->nsend; i++)
                        {
                            rvec_add(x[spas->a[i]], shift, *vbuf);
                            vbuf++;
                        }
                    }
                    else
                    {
                        /* Shift and rotate coordinates */
                        for (i = 0; i < spas->nsend; i++)
                        {
                            (*vbuf)[XX] =               x[spas->a[i]][XX] + shift[XX];
                            (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
                            (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
                            vbuf++;
                        }
                    }
                }
            }
            /* Send and receive the coordinates */
            spas = spac->spas[d];
            ns0  = spas[0].nsend;
            nr0  = spas[0].nrecv;
            ns1  = spas[1].nsend;
            nr1  = spas[1].nrecv;
            if (nvec == 1)
            {
                dd_sendrecv2_rvec(dd, d,
                                  spac->vbuf+ns0, ns1, x0+n, nr1,
                                  spac->vbuf, ns0, x0+n+nr1, nr0);
            }
            else
            {
                /* Communicate both vectors in one buffer */
                rbuf = spac->vbuf2;
                dd_sendrecv2_rvec(dd, d,
//.........这里部分代码省略.........