int n1, result, natoms;
Molecule *mol2;
IntGroup *ig;
- MolAction *act2;
+ MolAction *act2 = NULL;
int needsSymmetryAmendment = 0;
int needsRebuildMDArena = 0;
Int *ip;
} else if (strcmp(action->name, gMolActionAddBonds) == 0) {
ip = (Int *)action->args[0].u.arval.ptr;
n1 = action->args[0].u.arval.nitems / 2;
- if ((result = MoleculeAddBonds(mol, n1, ip)) < 0)
+ if ((result = MoleculeAddBonds(mol, n1, ip)) <= 0)
return result;
- act2 = MolActionNew(gMolActionDeleteBonds, n1 * 2, ip);
+ ip = (Int *)malloc(sizeof(Int) * 2 * result);
+ if (ip == NULL)
+ return -4;
+ memmove(ip, mol->bonds - result * 2, sizeof(Int) * 2 * result);
+ act2 = MolActionNew(gMolActionDeleteBonds, result * 2, ip);
+ free(ip);
needsRebuildMDArena = 1;
} else if (strcmp(action->name, gMolActionDeleteBonds) == 0) {
ip = (Int *)action->args[0].u.arval.ptr;
int
MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds)
{
- int i, j, n1, n2;
+ int i, j, n1, n2, n;
Atom *ap;
+ Int *bonds_tmp;
+
if (mp == NULL || bonds == NULL || nbonds <= 0)
return 0;
if (mp->noModifyTopology)
return -4; /* Prohibited operation */
/* Check the bonds */
+ bonds_tmp = (Int *)malloc(sizeof(Int) * nbonds * 2);
+ if (bonds_tmp == NULL)
+ return -4; /* Out of memory */
+ n = 0;
for (i = 0; i < nbonds; i++) {
n1 = bonds[i * 2];
n2 = bonds[i * 2 + 1];
ap = ATOM_AT_INDEX(mp->atoms, n1);
if (ap->nconnects >= ATOMS_MAX_CONNECTS - 1 || ATOM_AT_INDEX(mp->atoms, n2)->nconnects >= ATOMS_MAX_CONNECTS - 1)
return -2; /* Too many bonds */
+ /* Check duplicates */
for (j = 0; j < ap->nconnects; j++) {
if (ap->connects[j] == n2)
- return -3; /* Duplicate bond */
+ break;
}
+ if (j == ap->nconnects) {
+ bonds_tmp[n * 2] = n1;
+ bonds_tmp[n * 2 + 1] = n2;
+ n++;
+ }
+ }
+ if (n == 0) {
+ /* No bonds to add */
+ free(bonds_tmp);
+ return 0;
}
__MoleculeLock(mp);
/* Add connects[] */
- for (i = 0; i < nbonds; i++) {
- n1 = bonds[i * 2];
- n2 = bonds[i * 2 + 1];
+ for (i = 0; i < n; i++) {
+ n1 = bonds_tmp[i * 2];
+ n2 = bonds_tmp[i * 2 + 1];
ap = ATOM_AT_INDEX(mp->atoms, n1);
ap->connects[ap->nconnects++] = n2;
ap = ATOM_AT_INDEX(mp->atoms, n2);
n1 = mp->nbonds;
/* if (AssignArray(&(mp->bonds), &(mp->nbonds), sizeof(Int) * 2, mp->nbonds + nb - 1, NULL) == NULL
|| sInsertElementsToArrayAtPositions(mp->bonds, n1, bonds, nb, sizeof(Int) * 2, where) != 0) */
- if (AssignArray(&(mp->bonds), &(mp->nbonds), sizeof(Int) * 2, mp->nbonds + nbonds - 1, NULL) == NULL)
+ if (AssignArray(&(mp->bonds), &(mp->nbonds), sizeof(Int) * 2, mp->nbonds + n - 1, NULL) == NULL)
goto panic;
- memmove(mp->bonds + n1 * 2, bonds, sizeof(Int) * 2 * nbonds);
+ memmove(mp->bonds + n1 * 2, bonds_tmp, sizeof(Int) * 2 * n);
/* Add angles, dihedrals, impropers */
{
angles = dihedrals = impropers = NULL;
nangles = ndihedrals = nimpropers = 0;
- for (i = 0; i < nbonds; i++) {
- n1 = bonds[i * 2];
- n2 = bonds[i * 2 + 1];
+ for (i = 0; i < n; i++) {
+ n1 = bonds_tmp[i * 2];
+ n2 = bonds_tmp[i * 2 + 1];
ap1 = ATOM_AT_INDEX(mp->atoms, n1);
ap2 = ATOM_AT_INDEX(mp->atoms, n2);
/* Angles X-n1-n2 */
mp->needsMDRebuild = 1;
__MoleculeUnlock(mp);
- return nbonds;
+ free(bonds_tmp);
+ return n;
panic:
__MoleculeUnlock(mp);
*
* Find the angles from the bonds. Returns the number of angles newly created.
*/
+/*
static VALUE
s_Molecule_FindAngles(VALUE self)
{
}
return INT2NUM(nip);
}
-
+*/
/*
* call-seq:
* find_dihedrals -> Integer
*
* Find the dihedrals from the bonds. Returns the number of dihedrals newly created.
*/
+/*
static VALUE
s_Molecule_FindDihedrals(VALUE self)
{
}
return INT2NUM(nip);
}
+*/
/*
* call-seq:
/*
* call-seq:
- * create_bond(n1, n2, ...) -> Molecule
+ * create_bond(n1, n2, ...) -> Integer
*
- * Create bonds between atoms n1 and n2, n3 and n4, and so on. Returns self.
+ * Create bonds between atoms n1 and n2, n3 and n4, and so on. If the corresponding bond is already present for a particular pair,
+ * do nothing for that pair. Returns the number of bonds actually created.
* This operation is undoable.
*/
static VALUE
s_Molecule_CreateBond(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- Int i, *ip;
+ Int i, *ip, old_nbonds;
if (argc == 0)
rb_raise(rb_eMolbyError, "missing arguments");
if (argc % 2 != 0)
for (i = 0; i < argc; i++)
ip[i] = s_Molecule_AtomIndexFromValue(mol, argv[i]);
ip[argc] = kInvalidIndex;
-// i = MoleculeAddBonds(mol, ip, NULL);
+ old_nbonds = mol->nbonds;
i = MolActionCreateAndPerform(mol, gMolActionAddBonds, argc, ip);
if (i == -1)
rb_raise(rb_eMolbyError, "atom index out of range");
rb_raise(rb_eMolbyError, "duplicate bonds");
else if (i != 0)
rb_raise(rb_eMolbyError, "error in creating bonds");
- return self;
+ return INT2NUM(mol->nbonds - old_nbonds);
+}
+
+/*
+ * call-seq:
+ * molecule.remove_bonds(n1, n2, ...) -> Integer
+ *
+ * Remove bonds between atoms n1 and n2, n3 and n4, and so on. If the corresponding bond is not present for
+ * a particular pair, do nothing for that pair. Returns the number of bonds actually removed.
+ * This operation is undoable.
+ */
+static VALUE
+s_Molecule_RemoveBond(int argc, VALUE *argv, VALUE self)
+{
+ Molecule *mol;
+ Int i, *ip, old_nbonds;
+ if (argc == 0)
+ rb_raise(rb_eMolbyError, "missing arguments");
+ if (argc % 2 != 0)
+ rb_raise(rb_eMolbyError, "bonds should be specified by pairs of atom indices");
+ Data_Get_Struct(self, Molecule, mol);
+ ip = ALLOC_N(Int, argc + 1);
+ for (i = 0; i < argc; i++)
+ ip[i] = s_Molecule_AtomIndexFromValue(mol, argv[i]);
+ ip[argc] = kInvalidIndex;
+ old_nbonds = mol->nbonds;
+ MolActionCreateAndPerform(mol, gMolActionDeleteBonds, argc, ip);
+ return INT2NUM(old_nbonds - mol->nbonds);
}
/*
rb_define_method(rb_cMolecule, "ps_per_step", s_Molecule_PsPerStep, 0);
rb_define_method(rb_cMolecule, "ps_per_step=", s_Molecule_SetPsPerStep, 1);
- rb_define_method(rb_cMolecule, "find_angles", s_Molecule_FindAngles, 0);
- rb_define_method(rb_cMolecule, "find_dihedrals", s_Molecule_FindDihedrals, 0);
+/* rb_define_method(rb_cMolecule, "find_angles", s_Molecule_FindAngles, 0);
+ rb_define_method(rb_cMolecule, "find_dihedrals", s_Molecule_FindDihedrals, 0); */
rb_define_method(rb_cMolecule, "nresidues", s_Molecule_Nresidues, 0);
rb_define_method(rb_cMolecule, "nresidues=", s_Molecule_ChangeNresidues, 1);
rb_define_method(rb_cMolecule, "max_residue_number", s_Molecule_MaxResSeq, -1);
rb_define_method(rb_cMolecule, "create_atom", s_Molecule_CreateAnAtom, -1);
rb_define_method(rb_cMolecule, "duplicate_atom", s_Molecule_DuplicateAnAtom, -1);
rb_define_method(rb_cMolecule, "create_bond", s_Molecule_CreateBond, -1);
+ rb_define_alias(rb_cMolecule, "create_bonds", "create_bond");
+ rb_define_method(rb_cMolecule, "remove_bond", s_Molecule_RemoveBond, -1);
+ rb_define_alias(rb_cMolecule, "remove_bonds", "remove_bond");
rb_define_method(rb_cMolecule, "add_angle", s_Molecule_AddAngle, 3);
rb_define_method(rb_cMolecule, "remove_angle", s_Molecule_RemoveAngle, 3);
rb_define_method(rb_cMolecule, "add_dihedral", s_Molecule_AddDihedral, 4);