-/*
- * call-seq:
- * bonds_on_border(group = selection) -> Array of Array of two Integers
- *
- * Returns an array of bonds that connect an atom in the group and an atom out
- * of the group. The first atom in the bond always belongs to the group. If no
- * such bonds are present, an empty array is returned.
- */
-static VALUE
-s_Molecule_BondsOnBorder(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- IntGroup *ig, *bg;
- VALUE gval, retval;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &gval);
- if (gval == Qnil) {
- ig = MoleculeGetSelection(mol);
- if (ig != NULL)
- IntGroupRetain(ig);
- } else {
- ig = s_Molecule_AtomGroupFromValue(self, gval);
- }
- retval = rb_ary_new();
- if (ig == NULL)
- return retval;
- bg = MoleculeSearchBondsAcrossAtomGroup(mol, ig);
- if (bg != NULL) {
- IntGroupIterator iter;
- Int i;
- IntGroupIteratorInit(bg, &iter);
- while ((i = IntGroupIteratorNext(&iter)) >= 0) {
- /* The atoms at the border */
- Int n1, n2;
- n1 = mol->bonds[i * 2];
- n2 = mol->bonds[i * 2 + 1];
- if (IntGroupLookupPoint(ig, n1) < 0) {
- int w = n1;
- n1 = n2;
- n2 = w;
- if (IntGroupLookupPoint(ig, n1) < 0)
- continue; /* Actually this is an internal error */
- }
- rb_ary_push(retval, rb_ary_new3(2, INT2NUM(n1), INT2NUM(n2)));
- }
- IntGroupIteratorRelease(&iter);
- }
- IntGroupRelease(bg);
- IntGroupRelease(ig);
- return retval;
-}
-
-/*
- * call-seq:
- * translate(vec, group = nil) -> Molecule
- *
- * Translate the molecule by vec. If group is given, only atoms in the group are moved.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_Translate(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE vec, group;
- Vector v;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "11", &vec, &group);
- ig = (NIL_P(group) ? NULL : s_Molecule_AtomGroupFromValue(self, group));
- VectorFromValue(vec, &v);
-// MoleculeTranslate(mol, &v, ig);
- MolActionCreateAndPerform(mol, gMolActionTranslateAtoms, &v, ig);
- if (ig != NULL)
- IntGroupRelease(ig);
- return self;
-}
-
-/*
- * call-seq:
- * rotate(axis, angle, center = [0,0,0], group = nil) -> Molecule
- *
- * Rotate the molecule. The axis must not a zero vector. angle is given in degree.
- * If group is given, only atoms in the group are moved.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_Rotate(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- volatile VALUE aval, anval, cval, gval;
- Double angle;
- Vector av, cv;
- Transform tr;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "22", &aval, &anval, &cval, &gval);
- ig = (NIL_P(gval) ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
- angle = NUM2DBL(rb_Float(anval)) * kDeg2Rad;
- VectorFromValue(aval, &av);
- if (NIL_P(cval))
- cv.x = cv.y = cv.z = 0.0;
- else
- VectorFromValue(cval, &cv);
- if (TransformForRotation(tr, &av, angle, &cv))
- rb_raise(rb_eMolbyError, "rotation axis cannot be a zero vector");
- MolActionCreateAndPerform(mol, gMolActionTransformAtoms, &tr, ig);
- if (ig != NULL)
- IntGroupRelease(ig);
- return self;
-}
-
-/*
- * call-seq:
- * reflect(axis, center = [0,0,0], group = nil) -> Molecule
- *
- * Reflect the molecule by the plane which is perpendicular to axis and including center.
- * axis must not be a zero vector.
- * If group is given, only atoms in the group are moved.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_Reflect(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- volatile VALUE aval, cval, gval;
- Vector av, cv;
- Transform tr;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "12", &aval, &cval, &gval);
- ig = (NIL_P(gval) ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
- VectorFromValue(aval, &av);
- if (NIL_P(cval))
- cv.x = cv.y = cv.z = 0.0;
- else
- VectorFromValue(cval, &cv);
- if (TransformForReflection(tr, &av, &cv))
- rb_raise(rb_eMolbyError, "reflection axis cannot be a zero vector");
- MolActionCreateAndPerform(mol, gMolActionTransformAtoms, &tr, ig);
- if (ig != NULL)
- IntGroupRelease(ig);
- return self;
-}
-
-/*
- * call-seq:
- * invert(center = [0,0,0], group = nil) -> Molecule
- *
- * Invert the molecule with the given center.
- * If group is given, only atoms in the group are moved.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_Invert(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- volatile VALUE cval, gval;
- Vector cv;
- Transform tr;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "02", &cval, &gval);
- ig = (NIL_P(gval) ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
- if (NIL_P(cval))
- cv.x = cv.y = cv.z = 0.0;
- else
- VectorFromValue(cval, &cv);
- TransformForInversion(tr, &cv);
- MolActionCreateAndPerform(mol, gMolActionTransformAtoms, &tr, ig);
- if (ig != NULL)
- IntGroupRelease(ig);
- return self;
-}
-
-/*
- * call-seq:
- * transform(transform, group = nil) -> Molecule
- *
- * Transform the molecule by the given Transform object.
- * If group is given, only atoms in the group are moved.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_Transform(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE trans, group;
- Transform tr;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "11", &trans, &group);
- ig = (NIL_P(group) ? NULL : s_Molecule_AtomGroupFromValue(self, group));
- TransformFromValue(trans, &tr);
-/* MoleculeTransform(mol, tr, ig); */
- MolActionCreateAndPerform(mol, gMolActionTransformAtoms, &tr, ig);
- if (ig != NULL)
- IntGroupRelease(ig);
- return self;
-}
-
-static void
-s_Molecule_DoCenterOfMass(Molecule *mol, Vector *outv, IntGroup *ig)
-{
- switch (MoleculeCenterOfMass(mol, outv, ig)) {
- case 2: rb_raise(rb_eMolbyError, "atom group is empty"); break;
- case 3: rb_raise(rb_eMolbyError, "weight is zero --- atomic weights are not defined?"); break;
- case 0: break;
- default: rb_raise(rb_eMolbyError, "cannot calculate center of mass"); break;
- }
-}
-
-/*
- * call-seq:
- * center_of_mass(group = nil) -> Vector3D
- *
- * Calculate the center of mass for the given set of atoms. The argument
- * group is null, then all atoms are considered.
- */
-static VALUE
-s_Molecule_CenterOfMass(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE group;
- IntGroup *ig;
- Vector v;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &group);
- ig = (NIL_P(group) ? NULL : s_Molecule_AtomGroupFromValue(self, group));
- s_Molecule_DoCenterOfMass(mol, &v, ig);
- if (ig != NULL)
- IntGroupRelease(ig);
- return ValueFromVector(&v);
-}
-
-/*
- * call-seq:
- * centralize(group = nil) -> self
- *
- * Translate the molecule so that the center of mass of the given group is located
- * at (0, 0, 0). Equivalent to molecule.translate(molecule.center_of_mass(group) * -1).
- */
-static VALUE
-s_Molecule_Centralize(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE group;
- IntGroup *ig;
- Vector v;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &group);
- ig = (NIL_P(group) ? NULL : s_Molecule_AtomGroupFromValue(self, group));
- s_Molecule_DoCenterOfMass(mol, &v, ig);
- if (ig != NULL)
- IntGroupRelease(ig);
- v.x = -v.x;
- v.y = -v.y;
- v.z = -v.z;
- MolActionCreateAndPerform(mol, gMolActionTranslateAtoms, &v, NULL);
- return self;
-}
-
-/*
- * call-seq:
- * bounds(group = nil) -> [min, max]
- *
- * Calculate the boundary. The return value is an array of two Vector3D objects.
- */
-static VALUE
-s_Molecule_Bounds(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE group;
- IntGroup *ig;
- Vector vmin, vmax;
- int n;
- Atom *ap;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &group);
- ig = (NIL_P(group) ? NULL : s_Molecule_AtomGroupFromValue(self, group));
- if (ig != NULL && IntGroupGetCount(ig) == 0)
- rb_raise(rb_eMolbyError, "atom group is empty");
- vmin.x = vmin.y = vmin.z = 1e30;
- vmax.x = vmax.y = vmax.z = -1e30;
- for (n = 0, ap = mol->atoms; n < mol->natoms; n++, ap = ATOM_NEXT(ap)) {
- Vector r;
- if (ig != NULL && IntGroupLookup(ig, n, NULL) == 0)
- continue;
- r = ap->r;
- if (r.x < vmin.x)
- vmin.x = r.x;
- if (r.y < vmin.y)
- vmin.y = r.y;
- if (r.z < vmin.z)
- vmin.z = r.z;
- if (r.x > vmax.x)
- vmax.x = r.x;
- if (r.y > vmax.y)
- vmax.y = r.y;
- if (r.z > vmax.z)
- vmax.z = r.z;
- }
- return rb_ary_new3(2, ValueFromVector(&vmin), ValueFromVector(&vmax));
-}
-
-/* Get atom position or a vector */
-static void
-s_Molecule_GetVectorFromArg(Molecule *mol, VALUE val, Vector *vp)
-{
- if (rb_obj_is_kind_of(val, rb_cInteger) || rb_obj_is_kind_of(val, rb_cString)) {
- int n1 = s_Molecule_AtomIndexFromValue(mol, val);
- *vp = ATOM_AT_INDEX(mol->atoms, n1)->r;
- } else {
- VectorFromValue(val, vp);
- }
-}
-
-/*
- * call-seq:
- * measure_bond(n1, n2) -> Float
- *
- * Calculate the bond length. The arguments can either be atom indices, the "residue:name" representation,
- * or Vector3D values.
- * If the crystallographic cell is defined, then the internal coordinates are convereted to the cartesian.
- */
-static VALUE
-s_Molecule_MeasureBond(VALUE self, VALUE nval1, VALUE nval2)
-{
- Molecule *mol;
- Vector v1, v2;
- Data_Get_Struct(self, Molecule, mol);
- s_Molecule_GetVectorFromArg(mol, nval1, &v1);
- s_Molecule_GetVectorFromArg(mol, nval2, &v2);
- return rb_float_new(MoleculeMeasureBond(mol, &v1, &v2));
-}
-
-/*
- * call-seq:
- * measure_angle(n1, n2, n3) -> Float
- *
- * Calculate the bond angle. The arguments can either be atom indices, the "residue:name" representation,
- * or Vector3D values. The return value is in degree.
- * If the crystallographic cell is defined, then the internal coordinates are convereted to the cartesian.
- */
-static VALUE
-s_Molecule_MeasureAngle(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3)
-{
- Molecule *mol;
- Vector v1, v2, v3;
- Double d;
- Data_Get_Struct(self, Molecule, mol);
- s_Molecule_GetVectorFromArg(mol, nval1, &v1);
- s_Molecule_GetVectorFromArg(mol, nval2, &v2);
- s_Molecule_GetVectorFromArg(mol, nval3, &v3);
- d = MoleculeMeasureAngle(mol, &v1, &v2, &v3);
- if (isnan(d))
- return Qnil; /* Cannot define */
- else return rb_float_new(d);
-}
-
-/*
- * call-seq:
- * measure_dihedral(n1, n2, n3, n4) -> Float
- *
- * Calculate the dihedral angle. The arguments can either be atom indices, the "residue:name" representation,
- * or Vector3D values. The return value is in degree.
- * If the crystallographic cell is defined, then the internal coordinates are convereted to the cartesian.
- */
-static VALUE
-s_Molecule_MeasureDihedral(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3, VALUE nval4)
-{
- Molecule *mol;
- Vector v1, v2, v3, v4;
- Double d;
- Data_Get_Struct(self, Molecule, mol);
- s_Molecule_GetVectorFromArg(mol, nval1, &v1);
- s_Molecule_GetVectorFromArg(mol, nval2, &v2);
- s_Molecule_GetVectorFromArg(mol, nval3, &v3);
- s_Molecule_GetVectorFromArg(mol, nval4, &v4);
- d = MoleculeMeasureDihedral(mol, &v1, &v2, &v3, &v4);
- if (isnan(d))
- return Qnil; /* Cannot define */
- else return rb_float_new(d);
-}
-
-/*
- * call-seq:
- * expand_by_symmetry(group, sym, dx=0, dy=0, dz=0) -> Array
- *
- * Expand the specified part of the molecule by the given symmetry operation.
- * Returns the array of atom indices corresponding to the expanded atoms.
- */
-static VALUE
-s_Molecule_ExpandBySymmetry(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE gval, sval, xval, yval, zval, rval;
- IntGroup *ig;
- Int n[4];
- Int natoms;
- Int nidx, *idx;
-
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "23", &gval, &sval, &xval, &yval, &zval);
- n[0] = NUM2INT(rb_Integer(sval));
- n[1] = (xval == Qnil ? 0 : NUM2INT(rb_Integer(xval)));
- n[2] = (yval == Qnil ? 0 : NUM2INT(rb_Integer(yval)));
- n[3] = (zval == Qnil ? 0 : NUM2INT(rb_Integer(zval)));
- ig = s_Molecule_AtomGroupFromValue(self, gval);
- if (n[0] < 0 || (n[0] > 0 && n[0] >= mol->nsyms))
- rb_raise(rb_eMolbyError, "symmetry index is out of bounds");
- natoms = mol->natoms;
-
- MolActionCreateAndPerform(mol, gMolActionExpandBySymmetry, ig, n[1], n[2], n[3], n[0], &nidx, &idx);
-
- rval = rb_ary_new2(nidx);
- while (--nidx >= 0) {
- rb_ary_store(rval, nidx, INT2NUM(idx[nidx]));
- }
-/* if (natoms == mol->natoms)
- rval = Qnil;
- else {
- rval = IntGroup_Alloc(rb_cIntGroup);
- Data_Get_Struct(rval, IntGroup, ig);
- IntGroup_RaiseIfError(IntGroupAdd(ig, natoms, mol->natoms - natoms));
- } */
- return rval;
-}
-
-/*
- * call-seq:
- * amend_by_symmetry(group = nil) -> IntGroup
- *
- * Expand the specified part of the molecule by the given symmetry operation.
- * Returns an IntGroup containing the added atoms.
- */
-static VALUE
-s_Molecule_AmendBySymmetry(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- IntGroup *ig, *ig2;
- VALUE rval, gval;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &gval);
- if (gval != Qnil)
- ig = s_Molecule_AtomGroupFromValue(self, gval);
- else ig = NULL;
- MolActionCreateAndPerform(mol, gMolActionAmendBySymmetry, ig, &ig2);
- rval = ValueFromIntGroup(ig2);
- IntGroupRelease(ig2);
- return rval;
-}
-
-/*
- * call-seq:
- * transform_for_symop(symop, is_cartesian = nil) -> Transform
- *
- * Get the transform corresponding to the symmetry operation. The symop can either be
- * an integer (index of symmetry operation) or [sym, dx, dy, dz].
- * If is_cartesian is true, the returned transform is for cartesian coordinates.
- * Otherwise, the returned transform is for fractional coordinates.
- * Raises exception when no cell or no transform are defined.
- */
-static VALUE
-s_Molecule_TransformForSymop(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE sval, fval;
- Symop symop;
- Transform tr;
- Data_Get_Struct(self, Molecule, mol);
- if (mol->cell == NULL)
- rb_raise(rb_eMolbyError, "no unit cell is defined");
- if (mol->nsyms == 0)
- rb_raise(rb_eMolbyError, "no symmetry operation is defined");
- rb_scan_args(argc, argv, "11", &sval, &fval);
- if (rb_obj_is_kind_of(sval, rb_cNumeric)) {
- symop.sym = NUM2INT(rb_Integer(sval));
- symop.dx = symop.dy = symop.dz = 0;
- } else {
- sval = rb_ary_to_ary(sval);
- if (RARRAY_LEN(sval) < 4)
- rb_raise(rb_eMolbyError, "missing arguments as symop; at least four integers should be given");
- symop.sym = NUM2INT(rb_Integer(RARRAY_PTR(sval)[0]));
- symop.dx = NUM2INT(rb_Integer(RARRAY_PTR(sval)[1]));
- symop.dy = NUM2INT(rb_Integer(RARRAY_PTR(sval)[2]));
- symop.dz = NUM2INT(rb_Integer(RARRAY_PTR(sval)[3]));
- }
- if (symop.sym >= mol->nsyms)
- rb_raise(rb_eMolbyError, "index of symmetry operation (%d) is out of range", symop.sym);
- MoleculeGetTransformForSymop(mol, symop, &tr, (RTEST(fval) != 0));
- return ValueFromTransform(&tr);
-}
-
-/*
- * call-seq:
- * symop_for_transform(transform, is_cartesian = nil) -> [sym, dx, dy, dz]
- *
- * Get the symmetry operation corresponding to the given transform.
- * If is_cartesian is true, the given transform is for cartesian coordinates.
- * Otherwise, the given transform is for fractional coordinates.
- * Raises exception when no cell or no transform are defined.
- */
-static VALUE
-s_Molecule_SymopForTransform(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE tval, fval;
- Symop symop;
- Transform tr;
- int n;
- Data_Get_Struct(self, Molecule, mol);
- if (mol->cell == NULL)
- rb_raise(rb_eMolbyError, "no unit cell is defined");
- if (mol->nsyms == 0)
- rb_raise(rb_eMolbyError, "no symmetry operation is defined");
- rb_scan_args(argc, argv, "11", &tval, &fval);
- TransformFromValue(tval, &tr);
- n = MoleculeGetSymopForTransform(mol, tr, &symop, (RTEST(fval) != 0));
- if (n == 0) {
- return rb_ary_new3(4, INT2NUM(symop.sym), INT2NUM(symop.dx), INT2NUM(symop.dy), INT2NUM(symop.dz));
- } else {
- return Qnil; /* Not found */
- }
-}
-
-/*
- * call-seq:
- * wrap_unit_cell(group) -> Vector3D
- *
- * Move the specified group so that the center of mass of the group is within the
- * unit cell. The offset vector is returned. If no periodic box is defined,
- * exception is raised.
- */
-static VALUE
-s_Molecule_WrapUnitCell(VALUE self, VALUE gval)
-{
- Molecule *mol;
- IntGroup *ig;
- Vector v, cv, dv;
- Data_Get_Struct(self, Molecule, mol);
- if (mol->cell == NULL)
- rb_raise(rb_eMolbyError, "no unit cell is defined");
- ig = s_Molecule_AtomGroupFromValue(self, gval);
- s_Molecule_DoCenterOfMass(mol, &cv, ig);
- TransformVec(&v, mol->cell->rtr, &cv);
- if (mol->cell->flags[0])
- v.x -= floor(v.x);
- if (mol->cell->flags[1])
- v.y -= floor(v.y);
- if (mol->cell->flags[2])
- v.z -= floor(v.z);
- TransformVec(&dv, mol->cell->tr, &v);
- VecDec(dv, cv);
- MolActionCreateAndPerform(mol, gMolActionTranslateAtoms, &dv, ig);
- IntGroupRelease(ig);
- return ValueFromVector(&dv);
-}
-
-/*
- * call-seq:
- * find_conflicts(limit[, group1[, group2 [, ignore_exclusion]]]) -> [[n1, n2], [n3, n4], ...]
- *
- * Find pairs of atoms that are within the limit distance. If group1 and group2 are given, the
- * first and second atom in the pair should belong to group1 and group2, respectively.
- * If ignore_exclusion is true, then 1-2 (bonded), 1-3, 1-4 pairs are also considered.
- */
-static VALUE
-s_Molecule_FindConflicts(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE limval, gval1, gval2, rval, igval;
- IntGroup *ig1, *ig2;
- IntGroupIterator iter1, iter2;
- Int npairs, *pairs;
- Int n[2], i;
- Double lim;
- Vector r1;
- Atom *ap1, *ap2;
- MDExclusion *exinfo;
- Int *exlist;
-
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "13", &limval, &gval1, &gval2, &igval);
- lim = NUM2DBL(rb_Float(limval));
- if (lim <= 0.0)
- rb_raise(rb_eMolbyError, "the limit (%g) should be positive", lim);
- if (gval1 != Qnil)
- ig1 = s_Molecule_AtomGroupFromValue(self, gval1);
- else
- ig1 = IntGroupNewWithPoints(0, mol->natoms, -1);
- if (gval2 != Qnil)
- ig2 = s_Molecule_AtomGroupFromValue(self, gval2);
- else
- ig2 = IntGroupNewWithPoints(0, mol->natoms, -1);
-
- if (!RTEST(igval)) {
- /* Use the exclusion table in MDArena */
- if (mol->par == NULL || mol->arena == NULL || mol->arena->is_initialized == 0 || mol->needsMDRebuild) {
- VALUE mval = ValueFromMolecule(mol);
- s_RebuildMDParameterIfNecessary(mval, Qnil);
- }
- exinfo = mol->arena->exinfo; /* May be NULL */
- exlist = mol->arena->exlist;
- } else {
- exinfo = NULL;
- exlist = NULL;
- }
- IntGroupIteratorInit(ig1, &iter1);
- IntGroupIteratorInit(ig2, &iter2);
- npairs = 0;
- pairs = NULL;
- while ((n[0] = IntGroupIteratorNext(&iter1)) >= 0) {
- Int exn1, exn2;
- ap1 = ATOM_AT_INDEX(mol->atoms, n[0]);
- r1 = ap1->r;
- if (exinfo != NULL) {
- exn1 = exinfo[n[0]].index1;
- exn2 = exinfo[n[0] + 1].index1;
- } else exn1 = exn2 = -1;
- IntGroupIteratorReset(&iter2);
- while ((n[1] = IntGroupIteratorNext(&iter2)) >= 0) {
- ap2 = ATOM_AT_INDEX(mol->atoms, n[1]);
- if (n[0] == n[1])
- continue; /* Same atom */
- if (exinfo != NULL) {
- /* Look up exclusion table to exclude 1-2, 1-3, and 1-4 pairs */
- for (i = exn1; i < exn2; i++) {
- if (exlist[i] == n[1])
- break;
- }
- if (i < exn2)
- continue; /* Should be excluded */
- }
- if (MoleculeMeasureBond(mol, &r1, &(ap2->r)) < lim) {
- /* Is this pair already registered? */
- Int *ip;
- for (i = 0, ip = pairs; i < npairs; i++, ip += 2) {
- if ((ip[0] == n[0] && ip[1] == n[1]) || (ip[0] == n[1] && ip[1] == n[0]))
- break;
- }
- if (i >= npairs) {
- /* Not registered yet */
- AssignArray(&pairs, &npairs, sizeof(Int) * 2, npairs, n);
- }
- }
- }