#pragma mark ====== Molecule Class ======
+#pragma mark ------ Allocate/Release/Accessor ------
+
/* An malloc'ed string buffer. Retains the error/warning message from the last ***load/***save method. */
/* Accessible from Ruby as Molecule#error_message and Molecule#error_message=. */
char *gLoadSaveErrorMessage = NULL;
return ig;
}
-static void
-s_Molecule_RaiseOnLoadSave(int status, const char *msg, const char *fname)
-{
- if (gLoadSaveErrorMessage != NULL) {
- MyAppCallback_setConsoleColor(1);
- MyAppCallback_showScriptMessage("On loading %s:\n%s\n", fname, gLoadSaveErrorMessage);
- MyAppCallback_setConsoleColor(0);
- }
- if (status != 0)
- rb_raise(rb_eMolbyError, "%s %s", msg, fname);
-}
-
/*
* call-seq:
* dup -> Molecule
/*
* call-seq:
+ * atom_index(val) -> Integer
+ *
+ * Returns the atom index represented by val. val can be either a non-negative integer
+ * (directly representing the atom index), a negative integer (representing <code>natoms - val</code>),
+ * a string of type "resname.resid:name" or "resname:name" or "resid:name" or "name",
+ * where resname, resid, name are the residue name, residue id, and atom name respectively.
+ * If val is a string and multiple atoms match the description, the atom with the lowest index
+ * is returned.
+ */
+static VALUE
+s_Molecule_AtomIndex(VALUE self, VALUE val)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ return INT2NUM(s_Molecule_AtomIndexFromValue(mol, val));
+}
+
+/*
+ * call-seq:
+ * self == Molecule -> boolean
+ *
+ * True if the two arguments point to the same molecule.
+ */
+static VALUE
+s_Molecule_Equal(VALUE self, VALUE val)
+{
+ if (rb_obj_is_kind_of(val, rb_cMolecule)) {
+ Molecule *mol1, *mol2;
+ Data_Get_Struct(self, Molecule, mol1);
+ Data_Get_Struct(val, Molecule, mol2);
+ return (mol1 == mol2 ? Qtrue : Qfalse);
+ } else return Qfalse;
+}
+
+#pragma mark ------ Load/Save ------
+
+static void
+s_Molecule_RaiseOnLoadSave(int status, const char *msg, const char *fname)
+{
+ if (gLoadSaveErrorMessage != NULL) {
+ MyAppCallback_setConsoleColor(1);
+ MyAppCallback_showScriptMessage("On loading %s:\n%s\n", fname, gLoadSaveErrorMessage);
+ MyAppCallback_setConsoleColor(0);
+ }
+ if (status != 0)
+ rb_raise(rb_eMolbyError, "%s %s", msg, fname);
+}
+
+/*
+ * call-seq:
* loadmbsf(file) -> bool
*
* Read a structure from a mbsf file.
return Qtrue;
}
-/*
- * call-seq:
- * savetep(file) -> bool
- *
- * Write coordinates as an ORTEP file. Returns true if successful.
- */
-/*
-static VALUE
-s_Molecule_Savetep(VALUE self, VALUE fname)
-{
- char *fstr;
- Molecule *mol;
- char errbuf[128];
- Data_Get_Struct(self, Molecule, mol);
- fstr = FileStringValuePtr(fname);
- if (MoleculeWriteToTepFile(mol, fstr, errbuf, sizeof errbuf) != 0)
- rb_raise(rb_eMolbyError, errbuf);
- return Qtrue;
-}
-*/
-
/* load([ftype, ] fname, ...) */
static VALUE
s_Molecule_LoadSave(int argc, VALUE *argv, VALUE self, int loadFlag)
/*
* call-seq:
+ * open -> Molecule
+ * open(file) -> Molecule
+ *
+ * Create a new molecule from file as a document. If file is not given, an untitled document is created.
+ */
+static VALUE
+s_Molecule_Open(int argc, VALUE *argv, VALUE self)
+{
+ VALUE fname;
+ const char *p;
+ Molecule *mp;
+ VALUE iflag;
+ rb_scan_args(argc, argv, "01", &fname);
+ if (NIL_P(fname))
+ p = NULL;
+ else
+ p = FileStringValuePtr(fname);
+ iflag = Ruby_SetInterruptFlag(Qfalse);
+ mp = MoleculeCallback_openNewMolecule(p);
+ Ruby_SetInterruptFlag(iflag);
+ if (mp == NULL) {
+ if (p == NULL)
+ rb_raise(rb_eMolbyError, "Cannot create untitled document");
+ else
+ rb_raise(rb_eMolbyError, "Cannot open the file %s", p);
+ }
+ return ValueFromMolecule(mp);
+}
+
+/*
+ * call-seq:
+ * new -> Molecule
+ * new(file, *args) -> Molecule
+ *
+ * Create a new molecule and call "load" method with the same arguments.
+ */
+static VALUE
+s_Molecule_Initialize(int argc, VALUE *argv, VALUE self)
+{
+ if (argc > 0)
+ return s_Molecule_Load(argc, argv, self);
+ else return Qnil; /* An empty molecule (which is prepared in s_Molecule_Alloc()) is returned */
+}
+
+/*
+ * call-seq:
+ * error_message -> String
+ *
+ * Get the error_message from the last load/save method. If no error, returns nil.
+ */
+static VALUE
+s_Molecule_ErrorMessage(VALUE klass)
+{
+ if (gLoadSaveErrorMessage == NULL)
+ return Qnil;
+ else return rb_str_new2(gLoadSaveErrorMessage);
+}
+
+/*
+ * call-seq:
+ * set_error_message(String)
+ * Molecule.error_message = String
+ *
+ * Get the error_message from the last load/save method. If no error, returns nil.
+ */
+static VALUE
+s_Molecule_SetErrorMessage(VALUE klass, VALUE sval)
+{
+ if (gLoadSaveErrorMessage != NULL) {
+ free(gLoadSaveErrorMessage);
+ gLoadSaveErrorMessage = NULL;
+ }
+ if (sval != Qnil) {
+ sval = rb_str_to_str(sval);
+ gLoadSaveErrorMessage = strdup(StringValuePtr(sval));
+ }
+ return sval;
+}
+
+#pragma mark ------ Name attributes ------
+
+/*
+ * call-seq:
* name -> String
*
* Returns the display name of the molecule. If the molecule has no associated
}
}
-/*
- * call-seq:
- * open -> Molecule
- * open(file) -> Molecule
- *
- * Create a new molecule from file as a document. If file is not given, an untitled document is created.
- */
-static VALUE
-s_Molecule_Open(int argc, VALUE *argv, VALUE self)
-{
- VALUE fname;
- const char *p;
- Molecule *mp;
- VALUE iflag;
- rb_scan_args(argc, argv, "01", &fname);
- if (NIL_P(fname))
- p = NULL;
- else
- p = FileStringValuePtr(fname);
- iflag = Ruby_SetInterruptFlag(Qfalse);
- mp = MoleculeCallback_openNewMolecule(p);
- Ruby_SetInterruptFlag(iflag);
- if (mp == NULL) {
- if (p == NULL)
- rb_raise(rb_eMolbyError, "Cannot create untitled document");
- else
- rb_raise(rb_eMolbyError, "Cannot open the file %s", p);
- }
- return ValueFromMolecule(mp);
-}
-
-/*
- * call-seq:
- * new -> Molecule
- * new(file, *args) -> Molecule
- *
- * Create a new molecule and call "load" method with the same arguments.
- */
-static VALUE
-s_Molecule_Initialize(int argc, VALUE *argv, VALUE self)
-{
- if (argc > 0)
- return s_Molecule_Load(argc, argv, self);
- else return Qnil; /* An empty molecule (which is prepared in s_Molecule_Alloc()) is returned */
-}
+#pragma mark ------ MolEnumerables ------
static VALUE
s_Molecule_MolEnumerable(VALUE self, int kind)
return INT2NUM(mol->nresidues);
}
-static VALUE
-s_Molecule_BondParIsObsolete(VALUE self, VALUE val)
-{
- rb_raise(rb_eMolbyError, "Molecule#bond_par, angle_par, dihedral_par, improper_par, vdw_par are now obsolete. You can use MDArena#bond_par, angle_par, dihedral_par, improper_par, vdw_par instead, and probably these are what you really want.");
-}
-
/*
* call-seq:
- * bond_par(idx) -> ParameterRef
+ * nresidues = Integer
*
- * Returns the MD parameter for the idx-th bond.
+ * Change the number of residues.
*/
-/*
static VALUE
-s_Molecule_BondPar(VALUE self, VALUE val)
+s_Molecule_ChangeNresidues(VALUE self, VALUE val)
{
Molecule *mol;
- BondPar *bp;
- UInt t1, t2;
- Int i1, i2;
- Int ival;
- Data_Get_Struct(self, Molecule, mol);
- ival = NUM2INT(rb_Integer(val));
- if (ival < -mol->nbonds || ival >= mol->nbonds)
- rb_raise(rb_eMolbyError, "bond index (%d) out of range", ival);
- if (ival < 0)
- ival += mol->nbonds;
- s_RebuildMDParameterIfNecessary(self, Qtrue);
- i1 = mol->bonds[ival * 2];
- i2 = mol->bonds[ival * 2 + 1];
- t1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
- t2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
- bp = ParameterLookupBondPar(mol->par, t1, t2, i1, i2, 0);
- if (bp == NULL)
- return Qnil;
- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kBondParType, bp - mol->par->bondPars);
+ int ival = NUM2INT(val);
+ Data_Get_Struct(self, Molecule, mol);
+ MolActionCreateAndPerform(mol, gMolActionChangeNumberOfResidues, ival);
+ if (ival != mol->nresidues)
+ rb_raise(rb_eMolbyError, "Cannot set the number of residues to %d (set to %d)", ival, mol->nresidues);
+ return val;
}
-*/
/*
* call-seq:
- * angle_par(idx) -> ParameterRef
+ * max_residue_number(atom_group = nil) -> Integer
*
- * Returns the MD parameter for the idx-th angle.
+ * Returns the maximum residue number actually used. If an atom group is given, only
+ * these atoms are examined. If no atom is present, nil is returned.
*/
-/*
static VALUE
-s_Molecule_AnglePar(VALUE self, VALUE val)
+s_Molecule_MaxResSeq(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- AnglePar *ap;
- UInt t1, t2, t3;
- Int i1, i2, i3;
- Int ival;
- Data_Get_Struct(self, Molecule, mol);
- ival = NUM2INT(rb_Integer(val));
- if (ival < -mol->nangles || ival >= mol->nangles)
- rb_raise(rb_eMolbyError, "angle index (%d) out of range", ival);
- if (ival < 0)
- ival += mol->nangles;
- s_RebuildMDParameterIfNecessary(self, Qtrue);
- i1 = mol->angles[ival * 3];
- i2 = mol->angles[ival * 3 + 1];
- i3 = mol->angles[ival * 3 + 2];
- t1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
- t2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
- t3 = ATOM_AT_INDEX(mol->atoms, i3)->type;
- ap = ParameterLookupAnglePar(mol->par, t1, t2, t3, i1, i2, i3, 0);
- if (ap == NULL)
- return Qnil;
- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kAngleParType, ap - mol->par->anglePars);
+ VALUE gval;
+ int maxSeq;
+ IntGroup *ig;
+ Data_Get_Struct(self, Molecule, mol);
+ rb_scan_args(argc, argv, "01", &gval);
+ ig = (gval == Qnil ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
+ maxSeq = MoleculeMaximumResidueNumber(mol, ig);
+ return (maxSeq >= 0 ? INT2NUM(maxSeq) : Qnil);
}
-*/
+
/*
* call-seq:
- * dihedral_par(idx) -> ParameterRef
+ * min_residue_number(atom_group = nil) -> Integer
*
- * Returns the MD parameter for the idx-th dihedral.
+ * Returns the minimum residue number actually used. If an atom group is given, only
+ * these atoms are examined. If no atom is present, nil is returned.
*/
-/*
static VALUE
-s_Molecule_DihedralPar(VALUE self, VALUE val)
+s_Molecule_MinResSeq(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- Int ival;
- TorsionPar *tp;
- UInt t1, t2, t3, t4;
- Int i1, i2, i3, i4;
- Data_Get_Struct(self, Molecule, mol);
- ival = NUM2INT(rb_Integer(val));
- if (ival < -mol->ndihedrals || ival >= mol->ndihedrals)
- rb_raise(rb_eMolbyError, "dihedral index (%d) out of range", ival);
- if (ival < 0)
- ival += mol->ndihedrals;
- s_RebuildMDParameterIfNecessary(self, Qtrue);
- i1 = mol->dihedrals[ival * 4];
- i2 = mol->dihedrals[ival * 4 + 1];
- i3 = mol->dihedrals[ival * 4 + 2];
- i4 = mol->dihedrals[ival * 4 + 3];
- t1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
- t2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
- t3 = ATOM_AT_INDEX(mol->atoms, i3)->type;
- t4 = ATOM_AT_INDEX(mol->atoms, i4)->type;
- tp = ParameterLookupDihedralPar(mol->par, t1, t2, t3, t4, i1, i2, i3, i4, 0);
- if (tp == NULL)
- return Qnil;
- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kDihedralParType, tp - mol->par->dihedralPars);
+ VALUE gval;
+ int minSeq;
+ IntGroup *ig;
+ Data_Get_Struct(self, Molecule, mol);
+ rb_scan_args(argc, argv, "01", &gval);
+ ig = (gval == Qnil ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
+ minSeq = MoleculeMinimumResidueNumber(mol, ig);
+ return (minSeq >= 0 ? INT2NUM(minSeq) : Qnil);
}
-*/
+
/*
* call-seq:
- * improper_par(idx) -> ParameterRef
+ * each_atom(atom_group = nil) {|aref| ...}
*
- * Returns the MD parameter for the idx-th improper.
- */
-/*
-static VALUE
-s_Molecule_ImproperPar(VALUE self, VALUE val)
-{
- Molecule *mol;
- Int ival;
- TorsionPar *tp;
- UInt t1, t2, t3, t4;
- Int i1, i2, i3, i4;
- Data_Get_Struct(self, Molecule, mol);
- ival = NUM2INT(rb_Integer(val));
- if (ival < -mol->nimpropers || ival >= mol->nimpropers)
- rb_raise(rb_eMolbyError, "improper index (%d) out of range", ival);
- if (ival < 0)
- ival += mol->nimpropers;
- s_RebuildMDParameterIfNecessary(self, Qtrue);
- i1 = mol->impropers[ival * 4];
- i2 = mol->impropers[ival * 4 + 1];
- i3 = mol->impropers[ival * 4 + 2];
- i4 = mol->impropers[ival * 4 + 3];
- t1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
- t2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
- t3 = ATOM_AT_INDEX(mol->atoms, i3)->type;
- t4 = ATOM_AT_INDEX(mol->atoms, i4)->type;
- tp = ParameterLookupImproperPar(mol->par, t1, t2, t3, t4, i1, i2, i3, i4, 0);
- if (tp == NULL)
- return Qnil;
- return ValueFromMoleculeWithParameterTypeAndIndex(mol, kImproperParType, tp - mol->par->improperPars);
-}
-*/
-
-/*
- * call-seq:
- * start_step -> Integer
- *
- * Returns the start step (defined by dcd format).
+ * Execute the block, with the AtomRef object for each atom as the argument. If an atom
+ * group is given, only these atoms are processed.
+ * If atom_group is nil, this is equivalent to self.atoms.each, except that the return value
+ * is self (a Molecule object).
*/
static VALUE
-s_Molecule_StartStep(VALUE self)
+s_Molecule_EachAtom(int argc, VALUE *argv, VALUE self)
{
+ int i;
Molecule *mol;
+ AtomRef *aref;
+ VALUE arval;
+ VALUE gval;
+ IntGroup *ig;
Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(mol->startStep);
+ rb_scan_args(argc, argv, "01", &gval);
+ ig = (gval == Qnil ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
+ arval = ValueFromMoleculeAndIndex(mol, 0);
+ Data_Get_Struct(arval, AtomRef, aref);
+ for (i = 0; i < mol->natoms; i++) {
+ aref->idx = i;
+ if (ig == NULL || IntGroupLookup(ig, i, NULL))
+ rb_yield(arval);
+ }
+ if (ig != NULL)
+ IntGroupRelease(ig);
+ return self;
}
-/*
- * call-seq:
- * start_step = Integer
- *
- * Set the start step (defined by dcd format).
- */
+#pragma mark ------ Atom Group ------
+
static VALUE
-s_Molecule_SetStartStep(VALUE self, VALUE val)
+s_Molecule_AtomGroup_i(VALUE arg, VALUE values)
{
- Molecule *mol;
- Data_Get_Struct(self, Molecule, mol);
- mol->startStep = NUM2INT(rb_Integer(val));
- return val;
+ Molecule *mol = (Molecule *)(((VALUE *)values)[0]);
+ IntGroup *ig1 = (IntGroup *)(((VALUE *)values)[1]);
+ int idx = s_Molecule_AtomIndexFromValue(mol, arg);
+ IntGroup_RaiseIfError(IntGroupAdd(ig1, idx, 1));
+ return Qnil;
}
/*
* call-seq:
- * steps_per_frame -> Integer
+ * atom_group
+ * atom_group {|aref| ...}
+ * atom_group(arg1, arg2, ...)
+ * atom_group(arg1, arg2, ...) {|aref| ...}
+ *
+ * Specify a group of atoms. If no arguments are given, IntGroup\[0...natoms] is the result.
+ * If arguments are given, then the atoms reprensented by the arguments are added to the
+ * group. For a conversion of a string to an atom index, see the description
+ * of Molecule#atom_index.
+ * If a block is given, it is evaluated with an AtomRef (not atom index integers)
+ * representing each atom, and the atoms are removed from the result if the block returns false.
*
- * Returns the number of steps between frames (defined by dcd format).
*/
static VALUE
-s_Molecule_StepsPerFrame(VALUE self)
+s_Molecule_AtomGroup(int argc, VALUE *argv, VALUE self)
{
+ IntGroup *ig1, *ig2;
Molecule *mol;
+ Int i, startPt, interval;
+ VALUE retval = IntGroup_Alloc(rb_cIntGroup);
+ Data_Get_Struct(retval, IntGroup, ig1);
Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(mol->stepsPerFrame);
+ if (argc == 0) {
+ IntGroup_RaiseIfError(IntGroupAdd(ig1, 0, mol->natoms));
+ } else {
+ while (argc > 0) {
+ if (FIXNUM_P(*argv) || TYPE(*argv) == T_STRING) {
+ i = s_Molecule_AtomIndexFromValue(mol, *argv);
+ IntGroup_RaiseIfError(IntGroupAdd(ig1, i, 1));
+ } else if (rb_obj_is_kind_of(*argv, rb_cIntGroup)) {
+ ig2 = IntGroupFromValue(*argv);
+ for (i = 0; (startPt = IntGroupGetStartPoint(ig2, i)) >= 0; i++) {
+ interval = IntGroupGetInterval(ig2, i);
+ IntGroup_RaiseIfError(IntGroupAdd(ig1, startPt, interval));
+ }
+ IntGroupRelease(ig2);
+ } else if (rb_respond_to(*argv, rb_intern("each"))) {
+ VALUE values[2];
+ values[0] = (VALUE)mol;
+ values[1] = (VALUE)ig1;
+ rb_iterate(rb_each, *argv, s_Molecule_AtomGroup_i, (VALUE)values);
+ } else
+ IntGroup_RaiseIfError(IntGroupAdd(ig1, NUM2INT(*argv), 1));
+ argc--;
+ argv++;
+ }
+ }
+ if (rb_block_given_p()) {
+ /* Evaluate the given block with an AtomRef as the argument, and delete
+ the index if the block returns false */
+ AtomRef *aref = AtomRefNew(mol, 0);
+ VALUE arval = Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
+ ig2 = IntGroupNew();
+ IntGroupCopy(ig2, ig1);
+ for (i = 0; (startPt = IntGroupGetNthPoint(ig2, i)) >= 0; i++) {
+ VALUE resval;
+ if (startPt >= mol->natoms)
+ break;
+ aref->idx = startPt;
+ resval = rb_yield(arval);
+ if (!RTEST(resval))
+ IntGroupRemove(ig1, startPt, 1);
+ }
+ IntGroupRelease(ig2);
+ }
+
+ /* Remove points that are out of bounds */
+ IntGroup_RaiseIfError(IntGroupRemove(ig1, mol->natoms, INT_MAX));
+
+ return retval;
}
/*
* call-seq:
- * steps_per_frame = Integer
+ * selection -> IntGroup
*
- * Set the number of steps between frames (defined by dcd format).
+ * Returns the current selection.
*/
static VALUE
-s_Molecule_SetStepsPerFrame(VALUE self, VALUE val)
+s_Molecule_Selection(VALUE self)
{
Molecule *mol;
+ IntGroup *ig;
+ VALUE val;
Data_Get_Struct(self, Molecule, mol);
- mol->stepsPerFrame = NUM2INT(rb_Integer(val));
+ if (mol != NULL && (ig = MoleculeGetSelection(mol)) != NULL) {
+ ig = IntGroupNewFromIntGroup(ig); /* Duplicate, so that the change from GUI does not affect the value */
+ val = ValueFromIntGroup(ig);
+ IntGroupRelease(ig);
+ } else {
+ val = IntGroup_Alloc(rb_cIntGroup);
+ }
return val;
}
-/*
- * call-seq:
- * ps_per_step -> Float
- *
- * Returns the time increment (in picoseconds) for one step (defined by dcd format).
- */
static VALUE
-s_Molecule_PsPerStep(VALUE self)
+s_Molecule_SetSelectionSub(VALUE self, VALUE val, int undoable)
{
Molecule *mol;
+ IntGroup *ig;
Data_Get_Struct(self, Molecule, mol);
- return rb_float_new(mol->psPerStep);
+ if (val == Qnil)
+ ig = NULL;
+ else
+ ig = s_Molecule_AtomGroupFromValue(self, val);
+ if (undoable)
+ MolActionCreateAndPerform(mol, gMolActionSetSelection, ig);
+ else
+ MoleculeSetSelection(mol, ig);
+ if (ig != NULL)
+ IntGroupRelease(ig);
+ return val;
}
/*
* call-seq:
- * ps_per_step = Float
+ * selection = IntGroup
*
- * Set the time increment (in picoseconds) for one step (defined by dcd format).
+ * Set the current selection. The right-hand operand may be nil.
+ * This operation is _not_ undoable. If you need undo, use set_undoable_selection instead.
*/
static VALUE
-s_Molecule_SetPsPerStep(VALUE self, VALUE val)
+s_Molecule_SetSelection(VALUE self, VALUE val)
{
- Molecule *mol;
- Data_Get_Struct(self, Molecule, mol);
- mol->psPerStep = NUM2DBL(rb_Float(val));
- return val;
+ return s_Molecule_SetSelectionSub(self, val, 0);
}
/*
* call-seq:
- * find_angles -> Integer
+ * set_undoable_selection(IntGroup)
*
- * Find the angles from the bonds. Returns the number of angles newly created.
+ * Set the current selection with undo registration. The right-hand operand may be nil.
+ * This operation is undoable.
*/
-/*
static VALUE
-s_Molecule_FindAngles(VALUE self)
+s_Molecule_SetUndoableSelection(VALUE self, VALUE val)
{
- Molecule *mol;
- Atom *ap;
- int n1, i, j, nc;
- Int *ip, nip, n[3];
- Data_Get_Struct(self, Molecule, mol);
- if (mol == NULL || mol->natoms == 0)
- return INT2NUM(0);
- ip = NULL;
- nip = 0;
- for (n1 = 0, ap = mol->atoms; n1 < mol->natoms; n1++, ap = ATOM_NEXT(ap)) {
- nc = ap->connect.count;
- n[1] = n1;
- for (i = 0; i < nc; i++) {
- n[0] = ap->connects[i];
- for (j = i + 1; j < nc; j++) {
- n[2] = ap->connects[j];
- if (MoleculeLookupAngle(mol, n[0], n[1], n[2]) < 0)
- AssignArray(&ip, &nip, sizeof(Int) * 3, nip, n);
- }
- }
- }
- if (nip > 0) {
- MolActionCreateAndPerform(mol, gMolActionAddAngles, nip * 3, ip, NULL);
- free(ip);
- }
- return INT2NUM(nip);
+ return s_Molecule_SetSelectionSub(self, val, 1);
}
-*/
+
+#pragma mark ------ Editing ------
+
/*
* call-seq:
- * find_dihedrals -> Integer
+ * extract(group, dummy_flag = nil) -> Molecule
*
- * Find the dihedrals from the bonds. Returns the number of dihedrals newly created.
+ * Extract the atoms given by group and return as a new molecule object.
+ * If dummy_flag is true, then the atoms that are not included in the group but are connected
+ * to any atoms in the group are converted to "dummy" atoms (i.e. with element "Du" and
+ * names beginning with an underscore) and included in the new molecule object.
*/
-/*
static VALUE
-s_Molecule_FindDihedrals(VALUE self)
+s_Molecule_Extract(int argc, VALUE *argv, VALUE self)
{
- Molecule *mol;
- Atom *ap1, *ap2;
- int n1, i, j, k, nc1, nc2;
- Int *ip, nip, n[4];
- Data_Get_Struct(self, Molecule, mol);
- if (mol == NULL || mol->natoms == 0)
- return INT2NUM(0);
- ip = NULL;
- nip = 0;
- for (n1 = 0, ap1 = mol->atoms; n1 < mol->natoms; n1++, ap1 = ATOM_NEXT(ap1)) {
- nc1 = ap1->connect.count;
- n[1] = n1;
- for (i = 0; i < nc1; i++) {
- n[2] = ap1->connects[i];
- if (n[1] > n[2])
- continue;
- ap2 = ATOM_AT_INDEX(mol->atoms, n[2]);
- nc2 = ap2->connect.count;
- for (j = 0; j < nc1; j++) {
- n[0] = ap1->connects[j];
- if (n[0] == n[2])
- continue;
- for (k = 0; k < nc2; k++) {
- n[3] = ap2->connects[k];
- if (n[3] == n1 || n[3] == n[0])
- continue;
- if (MoleculeLookupDihedral(mol, n[0], n[1], n[2], n[3]) < 0)
- AssignArray(&ip, &nip, sizeof(Int) * 4, nip, n);
- }
- }
- }
- }
- if (nip > 0) {
- MolActionCreateAndPerform(mol, gMolActionAddDihedrals, nip * 4, ip, NULL);
- free(ip);
+ Molecule *mol1, *mol2;
+ IntGroup *ig;
+ VALUE group, dummy_flag, retval;
+ Data_Get_Struct(self, Molecule, mol1);
+ rb_scan_args(argc, argv, "11", &group, &dummy_flag);
+ ig = s_Molecule_AtomGroupFromValue(self, group);
+ if (MoleculeExtract(mol1, &mol2, ig, (dummy_flag != Qnil && dummy_flag != Qfalse)) != 0) {
+ retval = Qnil;
+ } else {
+ retval = ValueFromMolecule(mol2);
}
- return INT2NUM(nip);
+ IntGroupRelease(ig);
+ return retval;
}
-*/
/*
* call-seq:
- * nresidues = Integer
+ * add(molecule2) -> self
*
- * Change the number of residues.
+ * Combine two molecules. The residue numbers of the newly added atoms may be renumbered to avoid
+ conflicts.
+ This operation is undoable.
*/
static VALUE
-s_Molecule_ChangeNresidues(VALUE self, VALUE val)
+s_Molecule_Add(VALUE self, VALUE val)
{
- Molecule *mol;
- int ival = NUM2INT(val);
- Data_Get_Struct(self, Molecule, mol);
- MolActionCreateAndPerform(mol, gMolActionChangeNumberOfResidues, ival);
- if (ival != mol->nresidues)
- rb_raise(rb_eMolbyError, "Cannot set the number of residues to %d (set to %d)", ival, mol->nresidues);
- return val;
+ Molecule *mol1, *mol2;
+ Data_Get_Struct(self, Molecule, mol1);
+ mol2 = MoleculeFromValue(val);
+ MolActionCreateAndPerform(mol1, gMolActionMergeMolecule, mol2, NULL);
+ return self;
}
/*
* call-seq:
- * max_residue_number(atom_group = nil) -> Integer
+ * remove(group) -> Molecule
*
- * Returns the maximum residue number actually used. If an atom group is given, only
- * these atoms are examined. If no atom is present, nil is returned.
+ * The atoms designated by the given group are removed from the molecule.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_MaxResSeq(int argc, VALUE *argv, VALUE self)
+s_Molecule_Remove(VALUE self, VALUE group)
{
- Molecule *mol;
- VALUE gval;
- int maxSeq;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &gval);
- ig = (gval == Qnil ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
- maxSeq = MoleculeMaximumResidueNumber(mol, ig);
- return (maxSeq >= 0 ? INT2NUM(maxSeq) : Qnil);
-}
-
+ Molecule *mol1;
+ IntGroup *ig, *bg;
+ Int i;
+ IntGroupIterator iter;
+
+ Data_Get_Struct(self, Molecule, mol1);
+ group = rb_funcall(self, rb_intern("atom_group"), 1, group);
+ if (!rb_obj_is_kind_of(group, rb_cIntGroup))
+ rb_raise(rb_eMolbyError, "IntGroup instance is expected");
+ Data_Get_Struct(group, IntGroup, ig);
+
+ /* Remove the bonds between the two fragments */
+ /* (This is necessary for undo to work correctly) */
+ IntGroupIteratorInit(ig, &iter);
+ bg = NULL;
+ while ((i = IntGroupIteratorNext(&iter)) >= 0) {
+ Atom *ap = ATOM_AT_INDEX(mol1->atoms, i);
+ Int j, *cp;
+ cp = AtomConnectData(&ap->connect);
+ for (j = 0; j < ap->connect.count; j++) {
+ int n = cp[j];
+ if (!IntGroupLookup(ig, n, NULL)) {
+ /* bond i-n, i is in ig and n is not */
+ int k = MoleculeLookupBond(mol1, i, n);
+ if (k >= 0) {
+ if (bg == NULL)
+ bg = IntGroupNew();
+ IntGroupAdd(bg, k, 1);
+ }
+ }
+ }
+ }
+ IntGroupIteratorRelease(&iter);
+ if (bg != NULL) {
+ /* Remove bonds */
+ MolActionCreateAndPerform(mol1, gMolActionDeleteBonds, bg);
+ IntGroupRelease(bg);
+ }
+ /* Remove atoms */
+ if (MolActionCreateAndPerform(mol1, gMolActionUnmergeMolecule, ig) == 0)
+ return Qnil;
+ return self;
+}
+
/*
* call-seq:
- * min_residue_number(atom_group = nil) -> Integer
+ * create_atom(name, pos = -1) -> AtomRef
*
- * Returns the minimum residue number actually used. If an atom group is given, only
- * these atoms are examined. If no atom is present, nil is returned.
+ * Create a new atom with the specified name (may contain residue
+ * information) and position (if position is out of range, the atom is appended at
+ * the end). Returns the reference to the new atom.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_MinResSeq(int argc, VALUE *argv, VALUE self)
+s_Molecule_CreateAnAtom(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE gval;
- int minSeq;
- IntGroup *ig;
+ Int i, pos;
+ VALUE name, ival;
+ Atom arec;
+ AtomRef *aref;
+ char *p, resName[6], atomName[6];
+ int resSeq;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &gval);
- ig = (gval == Qnil ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
- minSeq = MoleculeMinimumResidueNumber(mol, ig);
- return (minSeq >= 0 ? INT2NUM(minSeq) : Qnil);
+ rb_scan_args(argc, argv, "02", &name, &ival);
+ if (ival != Qnil)
+ pos = NUM2INT(rb_Integer(ival));
+ else pos = -1;
+ if (name != Qnil) {
+ p = StringValuePtr(name);
+ if (p[0] != 0) {
+ i = MoleculeAnalyzeAtomName(p, resName, &resSeq, atomName);
+ if (atomName[0] == 0)
+ rb_raise(rb_eMolbyError, "bad atom name specification: %s", p);
+ }
+ } else p = NULL;
+ if (p == NULL || p[0] == 0) {
+ memset(atomName, 0, 4);
+ resSeq = -1;
+ }
+ memset(&arec, 0, sizeof(arec));
+ strncpy(arec.aname, atomName, 4);
+ if (resSeq >= 0) {
+ strncpy(arec.resName, resName, 4);
+ arec.resSeq = resSeq;
+ }
+ arec.occupancy = 1.0;
+// i = MoleculeCreateAnAtom(mol, &arec);
+ if (MolActionCreateAndPerform(mol, gMolActionAddAnAtom, &arec, pos, &pos) != 0)
+ return Qnil;
+ aref = AtomRefNew(mol, pos);
+ return Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
}
/*
* call-seq:
- * each_atom(atom_group = nil) {|aref| ...}
+ * duplicate_atom(atomref, pos = -1) -> AtomRef
*
- * Execute the block, with the AtomRef object for each atom as the argument. If an atom
- * group is given, only these atoms are processed.
- * If atom_group is nil, this is equivalent to self.atoms.each, except that the return value
- * is self (a Molecule object).
+ * Create a new atom with the same attributes (but no bonding information)
+ * with the specified atom. Returns the reference to the new atom.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_EachAtom(int argc, VALUE *argv, VALUE self)
+s_Molecule_DuplicateAnAtom(int argc, VALUE *argv, VALUE self)
{
- int i;
Molecule *mol;
+ const Atom *apsrc;
+ Atom arec;
AtomRef *aref;
- VALUE arval;
- VALUE gval;
- IntGroup *ig;
+ VALUE retval, aval, ival;
+ Int pos;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &gval);
- ig = (gval == Qnil ? NULL : s_Molecule_AtomGroupFromValue(self, gval));
- arval = ValueFromMoleculeAndIndex(mol, 0);
- Data_Get_Struct(arval, AtomRef, aref);
- for (i = 0; i < mol->natoms; i++) {
- aref->idx = i;
- if (ig == NULL || IntGroupLookup(ig, i, NULL))
- rb_yield(arval);
+ rb_scan_args(argc, argv, "11", &aval, &ival);
+ if (FIXNUM_P(aval)) {
+ int idx = NUM2INT(aval);
+ if (idx < 0 || idx >= mol->natoms)
+ rb_raise(rb_eMolbyError, "atom index out of range: %d", idx);
+ apsrc = ATOM_AT_INDEX(mol->atoms, idx);
+ } else {
+ apsrc = s_AtomFromValue(aval);
}
- if (ig != NULL)
- IntGroupRelease(ig);
- return self;
+ if (apsrc == NULL)
+ rb_raise(rb_eMolbyError, "bad atom specification");
+ if (ival != Qnil)
+ pos = NUM2INT(rb_Integer(ival));
+ else pos = -1;
+ AtomDuplicate(&arec, apsrc);
+ arec.connect.count = 0;
+ if (MolActionCreateAndPerform(mol, gMolActionAddAnAtom, &arec, pos, &pos) != 0)
+ retval = Qnil;
+ else {
+ aref = AtomRefNew(mol, pos);
+ retval = Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
+ }
+ AtomClean(&arec);
+ return retval;
}
/*
* call-seq:
- * cell -> [a, b, c, alpha, beta, gamma [, sig_a, sig_b, sig_c, sig_alpha, sig_beta, sig_gamma]]
+ * create_bond(n1, n2, ...) -> Integer
*
- * Returns the unit cell parameters. If cell is not set, returns nil.
+ * 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_Cell(VALUE self)
+s_Molecule_CreateBond(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- int i;
- VALUE val;
+ Int i, j, k, *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);
- if (mol->cell == NULL)
- return Qnil;
- val = rb_ary_new2(6);
- for (i = 0; i < 6; i++)
- rb_ary_push(val, rb_float_new(mol->cell->cell[i]));
- if (mol->cell->has_sigma) {
- for (i = 0; i < 6; i++) {
- rb_ary_push(val, rb_float_new(mol->cell->cellsigma[i]));
+ ip = ALLOC_N(Int, argc + 1);
+ for (i = j = 0; i < argc; i++, j++) {
+ ip[j] = s_Molecule_AtomIndexFromValue(mol, argv[i]);
+ if (i % 2 == 1) {
+ if (MoleculeLookupBond(mol, ip[j - 1], ip[j]) >= 0)
+ j -= 2; /* This bond is already present: skip it */
+ else {
+ for (k = 0; k < j - 1; k += 2) {
+ if ((ip[k] == ip[j - 1] && ip[k + 1] == ip[j]) || (ip[k + 1] == ip[j - 1] && ip[k] == ip[j])) {
+ j -= 2; /* The same entry is already in the argument */
+ break;
+ }
+ }
+ }
}
}
- return val;
+ old_nbonds = mol->nbonds;
+ if (j > 0) {
+ ip[j] = kInvalidIndex;
+ i = MolActionCreateAndPerform(mol, gMolActionAddBonds, j, ip, NULL);
+ } else i = 0;
+ xfree(ip);
+ if (i == -1)
+ rb_raise(rb_eMolbyError, "atom index out of range");
+ else if (i == -2)
+ rb_raise(rb_eMolbyError, "too many bonds");
+ else if (i == -3)
+ rb_raise(rb_eMolbyError, "duplicate bonds");
+ else if (i == -5)
+ rb_raise(rb_eMolbyError, "cannot create bond to itself");
+ else if (i != 0)
+ rb_raise(rb_eMolbyError, "error in creating bonds");
+ return INT2NUM(mol->nbonds - old_nbonds);
}
/*
* call-seq:
- * cell = [a, b, c, alpha, beta, gamma [, sig_a, sig_b, sig_c, sig_alpha, sig_beta, sig_gamma]]
- * set_cell([a, b, c, alpha, beta, gamma[, sig_a, sig_b, sig_c, sig_alpha, sig_beta, sig_gamma]], convert_coord = nil)
+ * molecule.remove_bonds(n1, n2, ...) -> Integer
*
- * Set the unit cell parameters. If the cell value is nil, then clear the current cell.
- If the given argument has 12 or more members, then the second half of the parameters represents the sigma values.
- This operation is undoable.
- Convert_coord is a flag to specify that the coordinates should be transformed so that the fractional coordinates remain the same.
+ * 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_SetCell(int argc, VALUE *argv, VALUE self)
+s_Molecule_RemoveBond(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE val, cval;
- int i, convert_coord, n;
- double d[12];
+ Int i, j, n[2];
+ IntGroup *bg;
+ 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);
- rb_scan_args(argc, argv, "11", &val, &cval);
- if (val == Qnil) {
- n = 0;
- } else {
- int len;
- val = rb_ary_to_ary(val);
- len = RARRAY_LEN(val);
- if (len >= 12) {
- n = 12;
- } else if (len >= 6) {
- n = 6;
- } else rb_raise(rb_eMolbyError, "too few members for cell parameters (6 or 12 required)");
- for (i = 0; i < n; i++)
- d[i] = NUM2DBL(rb_Float((RARRAY_PTR(val))[i]));
+ bg = NULL;
+ for (i = j = 0; i < argc; i++, j = 1 - j) {
+ n[j] = s_Molecule_AtomIndexFromValue(mol, argv[i]);
+ if (j == 1) {
+ Int k = MoleculeLookupBond(mol, n[0], n[1]);
+ if (k >= 0) {
+ if (bg == NULL)
+ bg = IntGroupNew();
+ IntGroupAdd(bg, k, 1);
+ }
+ }
}
- convert_coord = (RTEST(cval) ? 1 : 0);
- MolActionCreateAndPerform(mol, gMolActionSetCell, n, d, convert_coord);
- return val;
+ if (bg != NULL) {
+ MolActionCreateAndPerform(mol, gMolActionDeleteBonds, bg);
+ i = IntGroupGetCount(bg);
+ IntGroupRelease(bg);
+ } else i = 0;
+ return INT2NUM(i);
}
/*
* call-seq:
- * box -> [avec, bvec, cvec, origin, flags]
+ * assign_bond_order(idx, d1)
+ * assign_bond_orders(group, [d1, d2, ...])
*
- * Get the unit cell information in the form of a periodic bounding box.
- * Avec, bvec, cvec, origin are Vector3D objects, and flags is a 3-member array of
- * Integers which define whether the system is periodic along the axis.
- * If no unit cell is defined, nil is returned.
+ * Assign bond order. In the first form, the bond order of the idx-th bond is set to d1 (a Float value).
+ * In the second form, the bond orders at the indices in the group are set to d1, d2, etc.
+ * At present, the bond orders are only used in UFF parameter guess, and not in the MM/MD calculations.
+ * (This may change in the future)
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_Box(VALUE self)
+s_Molecule_AssignBondOrder(VALUE self, VALUE idxval, VALUE dval)
{
Molecule *mol;
- VALUE v[5], val;
+ IntGroup *ig;
Data_Get_Struct(self, Molecule, mol);
- if (mol == NULL || mol->cell == NULL)
- return Qnil;
- v[0] = ValueFromVector(&(mol->cell->axes[0]));
- v[1] = ValueFromVector(&(mol->cell->axes[1]));
- v[2] = ValueFromVector(&(mol->cell->axes[2]));
- v[3] = ValueFromVector(&(mol->cell->origin));
- v[4] = rb_ary_new3(3, INT2NUM(mol->cell->flags[0]), INT2NUM(mol->cell->flags[1]), INT2NUM(mol->cell->flags[2]));
- val = rb_ary_new4(5, v);
- return val;
+ if (rb_obj_is_kind_of(idxval, rb_cNumeric)) {
+ /* The first form */
+ Int idx = NUM2INT(rb_Integer(idxval));
+ Double d1 = NUM2DBL(rb_Float(dval));
+ if (idx < 0 || idx >= mol->nbonds)
+ rb_raise(rb_eMolbyError, "the bond index (%d) is out of bounds", idx);
+ ig = IntGroupNewWithPoints(idx, 1, -1);
+ MolActionCreateAndPerform(mol, gMolActionAssignBondOrders, 1, &d1, ig);
+ IntGroupRelease(ig);
+ } else {
+ Int i, n;
+ Double *dp;
+ ig = IntGroupFromValue(idxval);
+ n = IntGroupGetCount(ig);
+ if (n == 0)
+ rb_raise(rb_eMolbyError, "the bond index is empty");
+ dval = rb_ary_to_ary(dval);
+ dp = (Double *)calloc(sizeof(Double), n);
+ for (i = 0; i < RARRAY_LEN(dval) && i < n; i++) {
+ dp[i] = NUM2DBL(rb_Float(RARRAY_PTR(dval)[i]));
+ }
+ MolActionCreateAndPerform(mol, gMolActionAssignBondOrders, n, dp, ig);
+ free(dp);
+ IntGroupRelease(ig);
+ }
+ return self;
}
/*
* call-seq:
- * set_box(avec, bvec, cvec, origin = [0, 0, 0], flags = [1, 1, 1], convert_coordinates = nil)
- * set_box(d, origin = [0, 0, 0])
- * set_box
+ * get_bond_order(idx) -> Float
+ * get_bond_orders(group) -> Array
*
- * Set the unit cell parameters. Avec, bvec, and cvec can be either a Vector3D or a number.
- If it is a number, the x/y/z axis vector is multiplied with the given number and used
- as the box vector.
- Flags, if present, is a 3-member array of Integers defining whether the system is
- periodic along the axis.
- If convert_coordinates is true, then the coordinates are converted so that the fractional coordinates remain the same.
- In the second form, an isotropic box with cell-length d is set.
- In the third form, the existing box is cleared.
- Note: the sigma of the cell parameters is not cleared unless the periodic box itself is cleared.
+ * Get the bond order. In the first form, the bond order of the idx-th bond is returned.
+ * In the second form, the bond orders at the indices in the group are returned as an array.
+ * If no bond order information have been assigned, returns nil (the first form)
+ * or an empty array (the second form).
*/
static VALUE
-s_Molecule_SetBox(VALUE self, VALUE aval)
+s_Molecule_GetBondOrder(VALUE self, VALUE idxval)
{
Molecule *mol;
- VALUE v[6];
- static Vector ax[3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
- Vector vv[3];
- Vector origin = {0, 0, 0};
- char flags[3];
- Double d;
- int i, convertCoordinates = 0;
+ IntGroup *ig;
+ Double *dp;
+ VALUE retval;
+ Int i, n, numericArg;
Data_Get_Struct(self, Molecule, mol);
- if (aval == Qnil) {
- MolActionCreateAndPerform(mol, gMolActionClearBox);
- return self;
- }
- aval = rb_ary_to_ary(aval);
- for (i = 0; i < 6; i++) {
- if (i < RARRAY_LEN(aval))
- v[i] = (RARRAY_PTR(aval))[i];
- else v[i] = Qnil;
- }
- if (v[0] == Qnil) {
- MolActionCreateAndPerform(mol, gMolActionClearBox);
- return self;
- }
- if ((v[1] == Qnil || v[2] == Qnil) && rb_obj_is_kind_of(v[0], rb_cNumeric)) {
- d = NUM2DBL(rb_Float(v[0]));
- for (i = 0; i < 3; i++)
- VecScale(vv[i], ax[i], d);
- if (v[1] != Qnil)
- VectorFromValue(v[1], &origin);
- flags[0] = flags[1] = flags[2] = 1;
+ if (rb_obj_is_kind_of(idxval, rb_cNumeric)) {
+ /* The first form */
+ Int idx = NUM2INT(rb_Integer(idxval));
+ if (idx < 0 || idx >= mol->nbonds)
+ rb_raise(rb_eMolbyError, "the bond index (%d) is out of bounds", idx);
+ if (mol->bondOrders == NULL)
+ return Qnil;
+ ig = IntGroupNewWithPoints(idx, 1, -1);
+ n = 1;
+ numericArg = 1;
} else {
- for (i = 0; i < 3; i++) {
- if (v[i] == Qnil) {
- VecZero(vv[i]);
- } else if (rb_obj_is_kind_of(v[i], rb_cNumeric)) {
- d = NUM2DBL(rb_Float(v[i]));
- VecScale(vv[i], ax[i], d);
- } else {
- VectorFromValue(v[i], &vv[i]);
- }
- flags[i] = (VecLength2(vv[i]) > 0.0);
- }
- if (v[3] != Qnil)
- VectorFromValue(v[3], &origin);
- if (v[4] != Qnil) {
- for (i = 0; i < 3; i++) {
- VALUE val = Ruby_ObjectAtIndex(v[4], i);
- flags[i] = (NUM2INT(rb_Integer(val)) != 0);
- }
- }
- if (RTEST(v[5]))
- convertCoordinates = 1;
+ if (mol->bondOrders == NULL)
+ return rb_ary_new();
+ ig = IntGroupFromValue(idxval);
+ n = IntGroupGetCount(ig);
+ if (n == 0)
+ rb_raise(rb_eMolbyError, "the bond index is empty");
+ numericArg = 0;
}
- MolActionCreateAndPerform(mol, gMolActionSetBox, &(vv[0]), &(vv[1]), &(vv[2]), &origin, (flags[0] * 4 + flags[1] * 2 + flags[2]), convertCoordinates);
- return self;
+ dp = (Double *)calloc(sizeof(Double), n);
+ MoleculeGetBondOrders(mol, dp, ig);
+ if (numericArg)
+ retval = rb_float_new(dp[0]);
+ else {
+ retval = rb_ary_new();
+ for (i = 0; i < n; i++)
+ rb_ary_push(retval, rb_float_new(dp[i]));
+ }
+ free(dp);
+ IntGroupRelease(ig);
+ return retval;
}
/*
* call-seq:
- * cell_periodicity -> [n1, n2, n3]
+ * bond_exist?(idx1, idx2) -> bool
*
- * Get flags denoting whether the cell is periodic along the a/b/c axes. If the cell is not defined
- * nil is returned.
+ * Returns true if bond exists between atoms idx1 and idx2, otherwise returns false.
+ * Imaginary bonds between a pi-anchor and member atoms are not considered.
*/
static VALUE
-s_Molecule_CellPeriodicity(VALUE self)
+s_Molecule_BondExist(VALUE self, VALUE ival1, VALUE ival2)
{
- Molecule *mol;
+ Molecule *mol;
+ Int idx1, idx2, i;
+ Atom *ap;
+ Int *cp;
Data_Get_Struct(self, Molecule, mol);
- if (mol->cell == NULL)
- return Qnil;
- return rb_ary_new3(3, INT2FIX((int)mol->cell->flags[0]), INT2FIX((int)mol->cell->flags[1]), INT2FIX((int)mol->cell->flags[2]));
+ idx1 = NUM2INT(rb_Integer(ival1));
+ idx2 = NUM2INT(rb_Integer(ival2));
+ if (idx1 < 0 || idx1 >= mol->natoms || idx2 < 0 || idx2 >= mol->natoms)
+ rb_raise(rb_eMolbyError, "Atom index (%d or %d) out of range", idx1, idx2);
+ ap = ATOM_AT_INDEX(mol->atoms, idx1);
+ cp = AtomConnectData(&ap->connect);
+ for (i = 0; i < ap->connect.count; i++) {
+ if (cp[i] == idx2)
+ return Qtrue;
+ }
+ return Qfalse;
}
/*
* call-seq:
- * self.cell_periodicity = [n1, n2, n3] or Integer or nil
- * set_cell_periodicity = [n1, n2, n3] or Integer or nil
+ * add_angle(n1, n2, n3) -> Molecule
*
- * Set whether the cell is periodic along the a/b/c axes. If an integer is given as an argument,
- * its bits 2/1/0 (from the lowest) correspond to the a/b/c axes. Nil is equivalent to [0, 0, 0].
- * If cell is not defined, exception is raised.
+ * Add angle n1-n2-n3. Returns self. Usually, angles are automatically added
+ * when a bond is created, so it is rarely necessary to use this method explicitly.
* This operation is undoable.
*/
static VALUE
-s_Molecule_SetCellPeriodicity(VALUE self, VALUE arg)
+s_Molecule_AddAngle(VALUE self, VALUE v1, VALUE v2, VALUE v3)
{
+ Int n[4];
Molecule *mol;
- Int flag;
Data_Get_Struct(self, Molecule, mol);
- if (mol->cell == NULL)
- rb_raise(rb_eMolbyError, "periodic cell is not defined");
- if (arg == Qnil)
- flag = 0;
- else if (rb_obj_is_kind_of(arg, rb_cNumeric))
- flag = NUM2INT(rb_Integer(arg));
- else {
- Int i;
- VALUE arg0;
- arg = rb_ary_to_ary(arg);
- flag = 0;
- for (i = 0; i < 3 && i < RARRAY_LEN(arg); i++) {
- arg0 = RARRAY_PTR(arg)[i];
- if (arg0 != Qnil && arg0 != Qfalse && arg0 != INT2FIX(0))
- flag |= (1 << (2 - i));
- }
- }
- MolActionCreateAndPerform(mol, gMolActionSetCellPeriodicity, flag);
- return arg;
+ n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
+ n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
+ n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
+ if (MoleculeLookupAngle(mol, n[0], n[1], n[2]) >= 0)
+ rb_raise(rb_eMolbyError, "angle %d-%d-%d is already present", n[0], n[1], n[2]);
+ n[3] = kInvalidIndex;
+ MolActionCreateAndPerform(mol, gMolActionAddAngles, 3, n, NULL);
+ return self;
}
/*
* call-seq:
- * cell_flexibility -> bool
+ * remove_angle(n1, n2, n3) -> Molecule
*
- * Returns the unit cell is flexible or not
+ * Remove angle n1-n2-n3. Returns self. Usually, angles are automatically removed
+ * when a bond is removed, so it is rarely necessary to use this method explicitly.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_CellFlexibility(VALUE self)
+s_Molecule_RemoveAngle(VALUE self, VALUE v1, VALUE v2, VALUE v3)
{
- rb_warn("cell_flexibility is obsolete (unit cell is always frame dependent)");
- return Qtrue;
-/* Molecule *mol;
+ Int n[4];
+ Molecule *mol;
+ IntGroup *ig;
Data_Get_Struct(self, Molecule, mol);
- if (mol->cell == NULL)
- return Qfalse;
- if (mol->useFlexibleCell)
- return Qtrue;
- else return Qfalse; */
+ n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
+ n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
+ n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
+ if ((n[3] = MoleculeLookupAngle(mol, n[0], n[1], n[2])) < 0)
+ rb_raise(rb_eMolbyError, "angle %d-%d-%d is not present", n[0], n[1], n[2]);
+ ig = IntGroupNewWithPoints(n[3], 1, -1);
+ MolActionCreateAndPerform(mol, gMolActionDeleteAngles, ig);
+ IntGroupRelease(ig);
+ return self;
}
/*
* call-seq:
- * self.cell_flexibility = bool
- * set_cell_flexibility(bool)
+ * add_dihedral(n1, n2, n3, n4) -> Molecule
*
- * Change the unit cell is flexible or not
+ * Add dihedral n1-n2-n3-n4. Returns self. Usually, dihedrals are automatically added
+ * when a bond is created, so it is rarely necessary to use this method explicitly.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_SetCellFlexibility(VALUE self, VALUE arg)
+s_Molecule_AddDihedral(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
{
- rb_warn("set_cell_flexibility is obsolete (unit cell is always frame dependent)");
- return self;
-/* Molecule *mol;
+ Int n[5];
+ Molecule *mol;
Data_Get_Struct(self, Molecule, mol);
- MolActionCreateAndPerform(mol, gMolActionSetCellFlexibility, RTEST(arg) != 0);
- return self; */
+ n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
+ n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
+ n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
+ n[3] = s_Molecule_AtomIndexFromValue(mol, v4);
+ if (MoleculeLookupDihedral(mol, n[0], n[1], n[2], n[3]) >= 0)
+ rb_raise(rb_eMolbyError, "dihedral %d-%d-%d-%d is already present", n[0], n[1], n[2], n[3]);
+ n[4] = kInvalidIndex;
+ MolActionCreateAndPerform(mol, gMolActionAddDihedrals, 4, n, NULL);
+ return self;
}
/*
* call-seq:
- * cell_transform -> Transform
+ * remove_dihedral(n1, n2, n3, n4) -> Molecule
*
- * Get the transform matrix that converts internal coordinates to cartesian coordinates.
- * If cell is not defined, nil is returned.
+ * Remove dihedral n1-n2-n3-n4. Returns self. Usually, dihedrals are automatically removed
+ * when a bond is removed, so it is rarely necessary to use this method explicitly.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_CellTransform(VALUE self)
+s_Molecule_RemoveDihedral(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
{
+ Int n[5];
Molecule *mol;
+ IntGroup *ig;
Data_Get_Struct(self, Molecule, mol);
- if (mol == NULL || mol->cell == NULL)
- return Qnil;
- return ValueFromTransform(&(mol->cell->tr));
+ n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
+ n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
+ n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
+ n[3] = s_Molecule_AtomIndexFromValue(mol, v4);
+ if ((n[4] = MoleculeLookupDihedral(mol, n[0], n[1], n[2], n[3])) < 0)
+ rb_raise(rb_eMolbyError, "dihedral %d-%d-%d-%d is not present", n[0], n[1], n[2], n[3]);
+ ig = IntGroupNewWithPoints(n[4], 1, -1);
+ MolActionCreateAndPerform(mol, gMolActionDeleteDihedrals, ig);
+ IntGroupRelease(ig);
+ return self;
}
/*
* call-seq:
- * symmetry -> Array of Transforms
- * symmetries -> Array of Transforms
+ * add_improper(n1, n2, n3, n4) -> Molecule
*
- * Get the currently defined symmetry operations. If no symmetry operation is defined,
- * returns an empty array.
+ * Add dihedral n1-n2-n3-n4. Returns self. Unlike angles and dihedrals, impropers are
+ * not automatically added when a new bond is created, so this method is more useful than
+ * the angle/dihedral counterpart.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_Symmetry(VALUE self)
+s_Molecule_AddImproper(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
{
+ Int n[5];
Molecule *mol;
- VALUE val;
- int i;
Data_Get_Struct(self, Molecule, mol);
- if (mol->nsyms <= 0)
- return rb_ary_new();
- val = rb_ary_new2(mol->nsyms);
- for (i = 0; i < mol->nsyms; i++) {
- rb_ary_push(val, ValueFromTransform(&mol->syms[i]));
- }
- return val;
+ n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
+ n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
+ n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
+ n[3] = s_Molecule_AtomIndexFromValue(mol, v4);
+ if (MoleculeLookupImproper(mol, n[0], n[1], n[2], n[3]) >= 0)
+ rb_raise(rb_eMolbyError, "improper %d-%d-%d-%d is already present", n[0], n[1], n[2], n[3]);
+ n[4] = kInvalidIndex;
+ MolActionCreateAndPerform(mol, gMolActionAddImpropers, 4, n, NULL);
+ return self;
}
/*
* call-seq:
- * nsymmetries -> Integer
+ * remove_improper(n1, n2, n3, n4) -> Molecule
+ * remove_improper(intgroup) -> Molecule
*
- * Get the number of currently defined symmetry operations.
+ * Remove improper n1-n2-n3-n4, or the specified impropers (in indices) in IntGroup.
+ * Returns self. Unlike angles and dihedrals, impropers are
+ * not automatically added when a new bond is created, so this method is more useful than
+ * the angle/dihedral counterpart.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_Nsymmetries(VALUE self)
+s_Molecule_RemoveImproper(int argc, VALUE *argv, VALUE self)
{
+ Int n[5];
+ VALUE v1, v2, v3, v4;
Molecule *mol;
+ IntGroup *ig;
Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(mol->nsyms);
+ if (argc == 1) {
+ ig = IntGroupFromValue(argv[0]);
+ } else {
+ rb_scan_args(argc, argv, "40", &v1, &v2, &v3, &v4);
+ n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
+ n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
+ n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
+ n[3] = s_Molecule_AtomIndexFromValue(mol, v4);
+ if ((n[4] = MoleculeLookupImproper(mol, n[0], n[1], n[2], n[3])) < 0)
+ rb_raise(rb_eMolbyError, "improper %d-%d-%d-%d is not present", n[0], n[1], n[2], n[3]);
+ ig = IntGroupNewWithPoints(n[4], 1, -1);
+ }
+ MolActionCreateAndPerform(mol, gMolActionDeleteImpropers, ig);
+ IntGroupRelease(ig);
+ return self;
}
/*
* call-seq:
- * add_symmetry(Transform) -> Integer
+ * assign_residue(group, res) -> Molecule
*
- * Add a new symmetry operation. If no symmetry operation is defined and the
- * given argument is not an identity transform, then also add an identity
- * transform at the index 0.
- * Returns the total number of symmetries after operation.
+ * Assign the specified atoms as the given residue. res can either be an integer, "resname"
+ * or "resname.resno". When the residue number is not specified, the residue number of
+ * the first atom in the group is used.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_AddSymmetry(VALUE self, VALUE trans)
+s_Molecule_AssignResidue(VALUE self, VALUE range, VALUE res)
{
Molecule *mol;
- Transform tr;
+ IntGroup *ig;
+ char *p, *pp, buf[16];
+ Int resid, n;
+ Atom *ap;
Data_Get_Struct(self, Molecule, mol);
- TransformFromValue(trans, &tr);
- MolActionCreateAndPerform(mol, gMolActionAddSymmetryOperation, &tr);
- return INT2NUM(mol->nsyms);
+
+ /* Parse the argument res */
+ if (FIXNUM_P(res)) {
+ /* We can assume Fixnum here because Bignum is non-realistic as residue numbers */
+ resid = NUM2INT(res);
+ buf[0] = 0;
+ } else {
+ p = StringValuePtr(res);
+ pp = strchr(p, '.');
+ if (pp != NULL) {
+ resid = atoi(pp + 1);
+ n = pp - p;
+ } else {
+ resid = -1;
+ n = strlen(p);
+ }
+ if (n > sizeof buf - 1)
+ n = sizeof buf - 1;
+ strncpy(buf, p, n);
+ buf[n] = 0;
+ }
+ ig = s_Molecule_AtomGroupFromValue(self, range);
+ if (ig == NULL || IntGroupGetCount(ig) == 0)
+ return Qnil;
+
+ if (resid < 0) {
+ /* Use the residue number of the first specified atom */
+ n = IntGroupGetNthPoint(ig, 0);
+ if (n >= mol->natoms)
+ rb_raise(rb_eMolbyError, "Atom index (%d) out of range", n);
+ ap = ATOM_AT_INDEX(mol->atoms, n);
+ resid = ap->resSeq;
+ }
+ /* Change the residue number */
+ MolActionCreateAndPerform(mol, gMolActionChangeResidueNumber, ig, resid);
+ /* Change the residue name if necessary */
+ if (buf[0] != 0) {
+ /* Int seqs[2];
+ seqs[0] = resid;
+ seqs[1] = kInvalidIndex; */
+ MolActionCreateAndPerform(mol, gMolActionChangeResidueNames, 1, &resid, 4, buf);
+ }
+ IntGroupRelease(ig);
+ return self;
}
/*
* call-seq:
- * remove_symmetry(count = nil) -> Integer
- * remove_symmetries(count = nil) -> Integer
+ * offset_residue(group, offset) -> Molecule
*
- * Remove the specified number of symmetry operations. The last added ones are removed
- * first. If count is nil, then all symmetry operations are removed. Returns the
- * number of leftover symmetries.
+ * Offset the residue number of the specified atoms. If any of the residue number gets
+ * negative, then exception is thrown.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_RemoveSymmetry(int argc, VALUE *argv, VALUE self)
+s_Molecule_OffsetResidue(VALUE self, VALUE range, VALUE offset)
{
Molecule *mol;
- VALUE cval;
- int i, n;
+ IntGroup *ig;
+ int ofs, result;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &cval);
- if (cval == Qnil)
- n = mol->nsyms - 1;
- else {
- n = NUM2INT(rb_Integer(cval));
- if (n < 0 || n > mol->nsyms)
- rb_raise(rb_eMolbyError, "the given count of symops is out of range");
- if (n == mol->nsyms)
- n = mol->nsyms - 1;
- }
- for (i = 0; i < n; i++)
- MolActionCreateAndPerform(mol, gMolActionDeleteSymmetryOperation);
- return INT2NUM(mol->nsyms);
+ ig = s_Molecule_AtomGroupFromValue(self, range);
+ ofs = NUM2INT(offset);
+ result = MolActionCreateAndPerform(mol, gMolActionOffsetResidueNumbers, ig, ofs, -1);
+ if (result > 0)
+ rb_raise(rb_eMolbyError, "residue number of atom %d becomes negative", result - 1);
+ IntGroupRelease(ig);
+ return self;
}
+/*
+ * call-seq:
+ * renumber_atoms(array) -> IntGroup
+ *
+ * Change the order of atoms so that the atoms specified in the array argument appear
+ * in this order from the top of the molecule. The atoms that are not included in array
+ * are placed after these atoms, and these atoms are returned as an intGroup.
+ * This operation is undoable.
+ */
static VALUE
-s_Molecule_AtomGroup_i(VALUE arg, VALUE values)
+s_Molecule_RenumberAtoms(VALUE self, VALUE array)
{
- Molecule *mol = (Molecule *)(((VALUE *)values)[0]);
- IntGroup *ig1 = (IntGroup *)(((VALUE *)values)[1]);
- int idx = s_Molecule_AtomIndexFromValue(mol, arg);
- IntGroup_RaiseIfError(IntGroupAdd(ig1, idx, 1));
- return Qnil;
+ Molecule *mol;
+ Int *new2old;
+ IntGroup *ig;
+ int i, n;
+ VALUE *valp, retval;
+ Data_Get_Struct(self, Molecule, mol);
+ if (TYPE(array) != T_ARRAY)
+ array = rb_funcall(array, rb_intern("to_a"), 0);
+ n = RARRAY_LEN(array);
+ valp = RARRAY_PTR(array);
+ new2old = ALLOC_N(Int, n + 1);
+ for (i = 0; i < n; i++)
+ new2old[i] = s_Molecule_AtomIndexFromValue(mol, valp[i]);
+ new2old[i] = kInvalidIndex;
+ i = MolActionCreateAndPerform(mol, gMolActionRenumberAtoms, i, new2old);
+ if (i == 1)
+ rb_raise(rb_eMolbyError, "Atom index out of range");
+ else if (i == 2)
+ rb_raise(rb_eMolbyError, "Duplicate entry");
+ else if (i == 3)
+ rb_raise(rb_eMolbyError, "Internal inconsistency during atom renumbering");
+ retval = IntGroup_Alloc(rb_cIntGroup);
+ Data_Get_Struct(retval, IntGroup, ig);
+ if (mol->natoms > n)
+ IntGroup_RaiseIfError(IntGroupAdd(ig, n, mol->natoms - n));
+ xfree(new2old);
+ return retval;
}
/*
* call-seq:
- * atom_group
- * atom_group {|aref| ...}
- * atom_group(arg1, arg2, ...)
- * atom_group(arg1, arg2, ...) {|aref| ...}
- *
- * Specify a group of atoms. If no arguments are given, IntGroup\[0...natoms] is the result.
- * If arguments are given, then the atoms reprensented by the arguments are added to the
- * group. For a conversion of a string to an atom index, see the description
- * of Molecule#atom_index.
- * If a block is given, it is evaluated with an AtomRef (not atom index integers)
- * representing each atom, and the atoms are removed from the result if the block returns false.
+ * set_atom_attr(index, key, value)
*
+ * Set the atom attribute for the specified atom.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_AtomGroup(int argc, VALUE *argv, VALUE self)
+s_Molecule_SetAtomAttr(VALUE self, VALUE idx, VALUE key, VALUE val)
{
- IntGroup *ig1, *ig2;
- Molecule *mol;
- Int i, startPt, interval;
- VALUE retval = IntGroup_Alloc(rb_cIntGroup);
- Data_Get_Struct(retval, IntGroup, ig1);
+ Molecule *mol;
+ VALUE aref, oldval;
Data_Get_Struct(self, Molecule, mol);
- if (argc == 0) {
- IntGroup_RaiseIfError(IntGroupAdd(ig1, 0, mol->natoms));
- } else {
- while (argc > 0) {
- if (FIXNUM_P(*argv) || TYPE(*argv) == T_STRING) {
- i = s_Molecule_AtomIndexFromValue(mol, *argv);
- IntGroup_RaiseIfError(IntGroupAdd(ig1, i, 1));
- } else if (rb_obj_is_kind_of(*argv, rb_cIntGroup)) {
- ig2 = IntGroupFromValue(*argv);
- for (i = 0; (startPt = IntGroupGetStartPoint(ig2, i)) >= 0; i++) {
- interval = IntGroupGetInterval(ig2, i);
- IntGroup_RaiseIfError(IntGroupAdd(ig1, startPt, interval));
- }
- IntGroupRelease(ig2);
- } else if (rb_respond_to(*argv, rb_intern("each"))) {
- VALUE values[2];
- values[0] = (VALUE)mol;
- values[1] = (VALUE)ig1;
- rb_iterate(rb_each, *argv, s_Molecule_AtomGroup_i, (VALUE)values);
- } else
- IntGroup_RaiseIfError(IntGroupAdd(ig1, NUM2INT(*argv), 1));
- argc--;
- argv++;
- }
- }
- if (rb_block_given_p()) {
- /* Evaluate the given block with an AtomRef as the argument, and delete
- the index if the block returns false */
- AtomRef *aref = AtomRefNew(mol, 0);
- VALUE arval = Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
- ig2 = IntGroupNew();
- IntGroupCopy(ig2, ig1);
- for (i = 0; (startPt = IntGroupGetNthPoint(ig2, i)) >= 0; i++) {
- VALUE resval;
- if (startPt >= mol->natoms)
- break;
- aref->idx = startPt;
- resval = rb_yield(arval);
- if (!RTEST(resval))
- IntGroupRemove(ig1, startPt, 1);
- }
- IntGroupRelease(ig2);
- }
-
- /* Remove points that are out of bounds */
- IntGroup_RaiseIfError(IntGroupRemove(ig1, mol->natoms, INT_MAX));
-
- return retval;
+ aref = ValueFromMoleculeAndIndex(mol, s_Molecule_AtomIndexFromValue(mol, idx));
+ oldval = s_AtomRef_GetAttr(aref, key);
+ if (val == Qundef)
+ return oldval;
+ s_AtomRef_SetAttr(aref, key, val);
+ return val;
}
/*
* call-seq:
- * atom_index(val) -> Integer
+ * get_atom_attr(index, key)
*
- * Returns the atom index represented by val. val can be either a non-negative integer
- * (directly representing the atom index), a negative integer (representing <code>natoms - val</code>),
- * a string of type "resname.resid:name" or "resname:name" or "resid:name" or "name",
- * where resname, resid, name are the residue name, residue id, and atom name respectively.
- * If val is a string and multiple atoms match the description, the atom with the lowest index
- * is returned.
+ * Get the atom attribute for the specified atom.
*/
static VALUE
-s_Molecule_AtomIndex(VALUE self, VALUE val)
+s_Molecule_GetAtomAttr(VALUE self, VALUE idx, VALUE key)
{
- Molecule *mol;
- Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(s_Molecule_AtomIndexFromValue(mol, val));
+ return s_Molecule_SetAtomAttr(self, idx, key, Qundef);
}
+#pragma mark ------ Undo Support ------
+
/*
* call-seq:
- * extract(group, dummy_flag = nil) -> Molecule
+ * register_undo(script, *args)
*
- * Extract the atoms given by group and return as a new molecule object.
- * If dummy_flag is true, then the atoms that are not included in the group but are connected
- * to any atoms in the group are converted to "dummy" atoms (i.e. with element "Du" and
- * names beginning with an underscore) and included in the new molecule object.
+ * Register an undo operation with the current molecule.
*/
static VALUE
-s_Molecule_Extract(int argc, VALUE *argv, VALUE self)
+s_Molecule_RegisterUndo(int argc, VALUE *argv, VALUE self)
{
- Molecule *mol1, *mol2;
- IntGroup *ig;
- VALUE group, dummy_flag, retval;
- Data_Get_Struct(self, Molecule, mol1);
- rb_scan_args(argc, argv, "11", &group, &dummy_flag);
- ig = s_Molecule_AtomGroupFromValue(self, group);
- if (MoleculeExtract(mol1, &mol2, ig, (dummy_flag != Qnil && dummy_flag != Qfalse)) != 0) {
- retval = Qnil;
- } else {
- retval = ValueFromMolecule(mol2);
- }
- IntGroupRelease(ig);
- return retval;
+ Molecule *mol;
+ VALUE script, args;
+ MolAction *act;
+ Data_Get_Struct(self, Molecule, mol);
+ rb_scan_args(argc, argv, "1*", &script, &args);
+ act = MolActionNew(SCRIPT_ACTION("R"), StringValuePtr(script), args);
+ MolActionCallback_registerUndo(mol, act);
+ return script;
}
/*
* call-seq:
- * add(molecule2) -> self
+ * undo_enabled? -> bool
*
- * Combine two molecules. The residue numbers of the newly added atoms may be renumbered to avoid
- conflicts.
- This operation is undoable.
+ * Returns true if undo is enabled for this molecule; otherwise no.
*/
static VALUE
-s_Molecule_Add(VALUE self, VALUE val)
+s_Molecule_UndoEnabled(VALUE self)
{
- Molecule *mol1, *mol2;
- Data_Get_Struct(self, Molecule, mol1);
- mol2 = MoleculeFromValue(val);
- MolActionCreateAndPerform(mol1, gMolActionMergeMolecule, mol2, NULL);
- return self;
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ if (MolActionCallback_isUndoRegistrationEnabled(mol))
+ return Qtrue;
+ else return Qfalse;
}
/*
* call-seq:
- * remove(group) -> Molecule
+ * undo_enabled = bool
*
- * The atoms designated by the given group are removed from the molecule.
- * This operation is undoable.
+ * Enable or disable undo.
*/
static VALUE
-s_Molecule_Remove(VALUE self, VALUE group)
+s_Molecule_SetUndoEnabled(VALUE self, VALUE val)
{
- Molecule *mol1;
- IntGroup *ig, *bg;
- Int i;
- IntGroupIterator iter;
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ MolActionCallback_setUndoRegistrationEnabled(mol, (val != Qfalse && val != Qnil));
+ return val;
+}
- Data_Get_Struct(self, Molecule, mol1);
- group = rb_funcall(self, rb_intern("atom_group"), 1, group);
- if (!rb_obj_is_kind_of(group, rb_cIntGroup))
- rb_raise(rb_eMolbyError, "IntGroup instance is expected");
- Data_Get_Struct(group, IntGroup, ig);
+#pragma mark ------ Measure ------
- /* Remove the bonds between the two fragments */
- /* (This is necessary for undo to work correctly) */
- IntGroupIteratorInit(ig, &iter);
- bg = NULL;
- while ((i = IntGroupIteratorNext(&iter)) >= 0) {
- Atom *ap = ATOM_AT_INDEX(mol1->atoms, i);
- Int j, *cp;
- cp = AtomConnectData(&ap->connect);
- for (j = 0; j < ap->connect.count; j++) {
- int n = cp[j];
- if (!IntGroupLookup(ig, n, NULL)) {
- /* bond i-n, i is in ig and n is not */
- int k = MoleculeLookupBond(mol1, i, n);
- if (k >= 0) {
- if (bg == NULL)
- bg = IntGroupNew();
- IntGroupAdd(bg, k, 1);
- }
- }
- }
- }
- IntGroupIteratorRelease(&iter);
- if (bg != NULL) {
- /* Remove bonds */
- MolActionCreateAndPerform(mol1, gMolActionDeleteBonds, bg);
- IntGroupRelease(bg);
+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;
}
- /* Remove atoms */
- if (MolActionCreateAndPerform(mol1, gMolActionUnmergeMolecule, ig) == 0)
- return Qnil;
- return self;
}
/*
* call-seq:
- * create_atom(name, pos = -1) -> AtomRef
+ * center_of_mass(group = nil) -> Vector3D
*
- * Create a new atom with the specified name (may contain residue
- * information) and position (if position is out of range, the atom is appended at
- * the end). Returns the reference to the new atom.
- * This operation is undoable.
+ * 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_CreateAnAtom(int argc, VALUE *argv, VALUE self)
+s_Molecule_CenterOfMass(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- Int i, pos;
- VALUE name, ival;
- Atom arec;
- AtomRef *aref;
- char *p, resName[6], atomName[6];
- int resSeq;
+ VALUE group;
+ IntGroup *ig;
+ Vector v;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "02", &name, &ival);
- if (ival != Qnil)
- pos = NUM2INT(rb_Integer(ival));
- else pos = -1;
- if (name != Qnil) {
- p = StringValuePtr(name);
- if (p[0] != 0) {
- i = MoleculeAnalyzeAtomName(p, resName, &resSeq, atomName);
- if (atomName[0] == 0)
- rb_raise(rb_eMolbyError, "bad atom name specification: %s", p);
- }
- } else p = NULL;
- if (p == NULL || p[0] == 0) {
- memset(atomName, 0, 4);
- resSeq = -1;
- }
- memset(&arec, 0, sizeof(arec));
- strncpy(arec.aname, atomName, 4);
- if (resSeq >= 0) {
- strncpy(arec.resName, resName, 4);
- arec.resSeq = resSeq;
- }
- arec.occupancy = 1.0;
-// i = MoleculeCreateAnAtom(mol, &arec);
- if (MolActionCreateAndPerform(mol, gMolActionAddAnAtom, &arec, pos, &pos) != 0)
- return Qnil;
- aref = AtomRefNew(mol, pos);
- return Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
+ 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:
- * duplicate_atom(atomref, pos = -1) -> AtomRef
+ * centralize(group = nil) -> self
*
- * Create a new atom with the same attributes (but no bonding information)
- * with the specified atom. Returns the reference to the new atom.
- * This operation is undoable.
+ * 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_DuplicateAnAtom(int argc, VALUE *argv, VALUE self)
+s_Molecule_Centralize(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- const Atom *apsrc;
- Atom arec;
- AtomRef *aref;
- VALUE retval, aval, ival;
- Int pos;
+ VALUE group;
+ IntGroup *ig;
+ Vector v;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "11", &aval, &ival);
- if (FIXNUM_P(aval)) {
- int idx = NUM2INT(aval);
- if (idx < 0 || idx >= mol->natoms)
- rb_raise(rb_eMolbyError, "atom index out of range: %d", idx);
- apsrc = ATOM_AT_INDEX(mol->atoms, idx);
- } else {
- apsrc = s_AtomFromValue(aval);
- }
- if (apsrc == NULL)
- rb_raise(rb_eMolbyError, "bad atom specification");
- if (ival != Qnil)
- pos = NUM2INT(rb_Integer(ival));
- else pos = -1;
- AtomDuplicate(&arec, apsrc);
- arec.connect.count = 0;
- if (MolActionCreateAndPerform(mol, gMolActionAddAnAtom, &arec, pos, &pos) != 0)
- retval = Qnil;
- else {
- aref = AtomRefNew(mol, pos);
- retval = Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
- }
- AtomClean(&arec);
- return retval;
-}
-
-/*
- * call-seq:
- * create_bond(n1, n2, ...) -> Integer
- *
- * 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, j, k, *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 = j = 0; i < argc; i++, j++) {
- ip[j] = s_Molecule_AtomIndexFromValue(mol, argv[i]);
- if (i % 2 == 1) {
- if (MoleculeLookupBond(mol, ip[j - 1], ip[j]) >= 0)
- j -= 2; /* This bond is already present: skip it */
- else {
- for (k = 0; k < j - 1; k += 2) {
- if ((ip[k] == ip[j - 1] && ip[k + 1] == ip[j]) || (ip[k + 1] == ip[j - 1] && ip[k] == ip[j])) {
- j -= 2; /* The same entry is already in the argument */
- break;
- }
- }
- }
- }
- }
- old_nbonds = mol->nbonds;
- if (j > 0) {
- ip[j] = kInvalidIndex;
- i = MolActionCreateAndPerform(mol, gMolActionAddBonds, j, ip, NULL);
- } else i = 0;
- xfree(ip);
- if (i == -1)
- rb_raise(rb_eMolbyError, "atom index out of range");
- else if (i == -2)
- rb_raise(rb_eMolbyError, "too many bonds");
- else if (i == -3)
- rb_raise(rb_eMolbyError, "duplicate bonds");
- else if (i == -5)
- rb_raise(rb_eMolbyError, "cannot create bond to itself");
- else if (i != 0)
- rb_raise(rb_eMolbyError, "error in creating bonds");
- 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, j, n[2];
- IntGroup *bg;
- 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);
- bg = NULL;
- for (i = j = 0; i < argc; i++, j = 1 - j) {
- n[j] = s_Molecule_AtomIndexFromValue(mol, argv[i]);
- if (j == 1) {
- Int k = MoleculeLookupBond(mol, n[0], n[1]);
- if (k >= 0) {
- if (bg == NULL)
- bg = IntGroupNew();
- IntGroupAdd(bg, k, 1);
- }
- }
- }
- if (bg != NULL) {
- MolActionCreateAndPerform(mol, gMolActionDeleteBonds, bg);
- i = IntGroupGetCount(bg);
- IntGroupRelease(bg);
- } else i = 0;
- return INT2NUM(i);
-}
-
-/*
- * call-seq:
- * assign_bond_order(idx, d1)
- * assign_bond_orders(group, [d1, d2, ...])
- *
- * Assign bond order. In the first form, the bond order of the idx-th bond is set to d1 (a Float value).
- * In the second form, the bond orders at the indices in the group are set to d1, d2, etc.
- * At present, the bond orders are only used in UFF parameter guess, and not in the MM/MD calculations.
- * (This may change in the future)
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_AssignBondOrder(VALUE self, VALUE idxval, VALUE dval)
-{
- Molecule *mol;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- if (rb_obj_is_kind_of(idxval, rb_cNumeric)) {
- /* The first form */
- Int idx = NUM2INT(rb_Integer(idxval));
- Double d1 = NUM2DBL(rb_Float(dval));
- if (idx < 0 || idx >= mol->nbonds)
- rb_raise(rb_eMolbyError, "the bond index (%d) is out of bounds", idx);
- ig = IntGroupNewWithPoints(idx, 1, -1);
- MolActionCreateAndPerform(mol, gMolActionAssignBondOrders, 1, &d1, ig);
- IntGroupRelease(ig);
- } else {
- Int i, n;
- Double *dp;
- ig = IntGroupFromValue(idxval);
- n = IntGroupGetCount(ig);
- if (n == 0)
- rb_raise(rb_eMolbyError, "the bond index is empty");
- dval = rb_ary_to_ary(dval);
- dp = (Double *)calloc(sizeof(Double), n);
- for (i = 0; i < RARRAY_LEN(dval) && i < n; i++) {
- dp[i] = NUM2DBL(rb_Float(RARRAY_PTR(dval)[i]));
- }
- MolActionCreateAndPerform(mol, gMolActionAssignBondOrders, n, dp, ig);
- free(dp);
+ 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:
- * get_bond_order(idx) -> Float
- * get_bond_orders(group) -> Array
+ * bounds(group = nil) -> [min, max]
*
- * Get the bond order. In the first form, the bond order of the idx-th bond is returned.
- * In the second form, the bond orders at the indices in the group are returned as an array.
- * If no bond order information have been assigned, returns nil (the first form)
- * or an empty array (the second form).
+ * Calculate the boundary. The return value is an array of two Vector3D objects.
*/
static VALUE
-s_Molecule_GetBondOrder(VALUE self, VALUE idxval)
+s_Molecule_Bounds(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
+ VALUE group;
IntGroup *ig;
- Double *dp;
- VALUE retval;
- Int i, n, numericArg;
- Data_Get_Struct(self, Molecule, mol);
- if (rb_obj_is_kind_of(idxval, rb_cNumeric)) {
- /* The first form */
- Int idx = NUM2INT(rb_Integer(idxval));
- if (idx < 0 || idx >= mol->nbonds)
- rb_raise(rb_eMolbyError, "the bond index (%d) is out of bounds", idx);
- if (mol->bondOrders == NULL)
- return Qnil;
- ig = IntGroupNewWithPoints(idx, 1, -1);
- n = 1;
- numericArg = 1;
- } else {
- if (mol->bondOrders == NULL)
- return rb_ary_new();
- ig = IntGroupFromValue(idxval);
- n = IntGroupGetCount(ig);
- if (n == 0)
- rb_raise(rb_eMolbyError, "the bond index is empty");
- numericArg = 0;
- }
- dp = (Double *)calloc(sizeof(Double), n);
- MoleculeGetBondOrders(mol, dp, ig);
- if (numericArg)
- retval = rb_float_new(dp[0]);
- else {
- retval = rb_ary_new();
- for (i = 0; i < n; i++)
- rb_ary_push(retval, rb_float_new(dp[i]));
- }
- free(dp);
- IntGroupRelease(ig);
- return retval;
-}
-
-/*
- * call-seq:
- * bond_exist?(idx1, idx2) -> bool
- *
- * Returns true if bond exists between atoms idx1 and idx2, otherwise returns false.
- * Imaginary bonds between a pi-anchor and member atoms are not considered.
- */
-static VALUE
-s_Molecule_BondExist(VALUE self, VALUE ival1, VALUE ival2)
-{
- Molecule *mol;
- Int idx1, idx2, i;
+ Vector vmin, vmax;
+ int n;
Atom *ap;
- Int *cp;
- Data_Get_Struct(self, Molecule, mol);
- idx1 = NUM2INT(rb_Integer(ival1));
- idx2 = NUM2INT(rb_Integer(ival2));
- if (idx1 < 0 || idx1 >= mol->natoms || idx2 < 0 || idx2 >= mol->natoms)
- rb_raise(rb_eMolbyError, "Atom index (%d or %d) out of range", idx1, idx2);
- ap = ATOM_AT_INDEX(mol->atoms, idx1);
- cp = AtomConnectData(&ap->connect);
- for (i = 0; i < ap->connect.count; i++) {
- if (cp[i] == idx2)
- return Qtrue;
- }
- return Qfalse;
-}
-
-/*
- * call-seq:
- * add_angle(n1, n2, n3) -> Molecule
- *
- * Add angle n1-n2-n3. Returns self. Usually, angles are automatically added
- * when a bond is created, so it is rarely necessary to use this method explicitly.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_AddAngle(VALUE self, VALUE v1, VALUE v2, VALUE v3)
-{
- Int n[4];
- Molecule *mol;
Data_Get_Struct(self, Molecule, mol);
- n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
- n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
- n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
- if (MoleculeLookupAngle(mol, n[0], n[1], n[2]) >= 0)
- rb_raise(rb_eMolbyError, "angle %d-%d-%d is already present", n[0], n[1], n[2]);
- n[3] = kInvalidIndex;
- MolActionCreateAndPerform(mol, gMolActionAddAngles, 3, n, NULL);
- return self;
-}
-
-/*
- * call-seq:
- * remove_angle(n1, n2, n3) -> Molecule
- *
- * Remove angle n1-n2-n3. Returns self. Usually, angles are automatically removed
- * when a bond is removed, so it is rarely necessary to use this method explicitly.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_RemoveAngle(VALUE self, VALUE v1, VALUE v2, VALUE v3)
-{
- Int n[4];
- Molecule *mol;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
- n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
- n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
- if ((n[3] = MoleculeLookupAngle(mol, n[0], n[1], n[2])) < 0)
- rb_raise(rb_eMolbyError, "angle %d-%d-%d is not present", n[0], n[1], n[2]);
- ig = IntGroupNewWithPoints(n[3], 1, -1);
- MolActionCreateAndPerform(mol, gMolActionDeleteAngles, ig);
- IntGroupRelease(ig);
- return self;
+ 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));
}
-/*
- * call-seq:
- * add_dihedral(n1, n2, n3, n4) -> Molecule
- *
- * Add dihedral n1-n2-n3-n4. Returns self. Usually, dihedrals are automatically added
- * when a bond is created, so it is rarely necessary to use this method explicitly.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_AddDihedral(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
+/* Get atom position or a vector */
+static void
+s_Molecule_GetVectorFromArg(Molecule *mol, VALUE val, Vector *vp)
{
- Int n[5];
- Molecule *mol;
- Data_Get_Struct(self, Molecule, mol);
- n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
- n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
- n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
- n[3] = s_Molecule_AtomIndexFromValue(mol, v4);
- if (MoleculeLookupDihedral(mol, n[0], n[1], n[2], n[3]) >= 0)
- rb_raise(rb_eMolbyError, "dihedral %d-%d-%d-%d is already present", n[0], n[1], n[2], n[3]);
- n[4] = kInvalidIndex;
- MolActionCreateAndPerform(mol, gMolActionAddDihedrals, 4, n, NULL);
- return self;
+ 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:
- * remove_dihedral(n1, n2, n3, n4) -> Molecule
+ * measure_bond(n1, n2) -> Float
*
- * Remove dihedral n1-n2-n3-n4. Returns self. Usually, dihedrals are automatically removed
- * when a bond is removed, so it is rarely necessary to use this method explicitly.
- * This operation is undoable.
+ * 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_RemoveDihedral(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
+s_Molecule_MeasureBond(VALUE self, VALUE nval1, VALUE nval2)
{
- Int n[5];
Molecule *mol;
- IntGroup *ig;
+ Vector v1, v2;
Data_Get_Struct(self, Molecule, mol);
- n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
- n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
- n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
- n[3] = s_Molecule_AtomIndexFromValue(mol, v4);
- if ((n[4] = MoleculeLookupDihedral(mol, n[0], n[1], n[2], n[3])) < 0)
- rb_raise(rb_eMolbyError, "dihedral %d-%d-%d-%d is not present", n[0], n[1], n[2], n[3]);
- ig = IntGroupNewWithPoints(n[4], 1, -1);
- MolActionCreateAndPerform(mol, gMolActionDeleteDihedrals, ig);
- IntGroupRelease(ig);
- return self;
+ s_Molecule_GetVectorFromArg(mol, nval1, &v1);
+ s_Molecule_GetVectorFromArg(mol, nval2, &v2);
+ return rb_float_new(MoleculeMeasureBond(mol, &v1, &v2));
}
/*
* call-seq:
- * add_improper(n1, n2, n3, n4) -> Molecule
+ * measure_angle(n1, n2, n3) -> Float
*
- * Add dihedral n1-n2-n3-n4. Returns self. Unlike angles and dihedrals, impropers are
- * not automatically added when a new bond is created, so this method is more useful than
- * the angle/dihedral counterpart.
- * This operation is undoable.
+ * 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_AddImproper(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
+s_Molecule_MeasureAngle(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3)
{
- Int n[5];
Molecule *mol;
+ Vector v1, v2, v3;
+ Double d;
Data_Get_Struct(self, Molecule, mol);
- n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
- n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
- n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
- n[3] = s_Molecule_AtomIndexFromValue(mol, v4);
- if (MoleculeLookupImproper(mol, n[0], n[1], n[2], n[3]) >= 0)
- rb_raise(rb_eMolbyError, "improper %d-%d-%d-%d is already present", n[0], n[1], n[2], n[3]);
- n[4] = kInvalidIndex;
- MolActionCreateAndPerform(mol, gMolActionAddImpropers, 4, n, NULL);
- return self;
+ 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:
- * remove_improper(n1, n2, n3, n4) -> Molecule
- * remove_improper(intgroup) -> Molecule
+ * measure_dihedral(n1, n2, n3, n4) -> Float
*
- * Remove improper n1-n2-n3-n4, or the specified impropers (in indices) in IntGroup.
- * Returns self. Unlike angles and dihedrals, impropers are
- * not automatically added when a new bond is created, so this method is more useful than
- * the angle/dihedral counterpart.
- * This operation is undoable.
+ * 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_RemoveImproper(int argc, VALUE *argv, VALUE self)
+s_Molecule_MeasureDihedral(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3, VALUE nval4)
{
- Int n[5];
- VALUE v1, v2, v3, v4;
Molecule *mol;
- IntGroup *ig;
+ Vector v1, v2, v3, v4;
+ Double d;
Data_Get_Struct(self, Molecule, mol);
- if (argc == 1) {
- ig = IntGroupFromValue(argv[0]);
- } else {
- rb_scan_args(argc, argv, "40", &v1, &v2, &v3, &v4);
- n[0] = s_Molecule_AtomIndexFromValue(mol, v1);
- n[1] = s_Molecule_AtomIndexFromValue(mol, v2);
- n[2] = s_Molecule_AtomIndexFromValue(mol, v3);
- n[3] = s_Molecule_AtomIndexFromValue(mol, v4);
- if ((n[4] = MoleculeLookupImproper(mol, n[0], n[1], n[2], n[3])) < 0)
- rb_raise(rb_eMolbyError, "improper %d-%d-%d-%d is not present", n[0], n[1], n[2], n[3]);
- ig = IntGroupNewWithPoints(n[4], 1, -1);
- }
- MolActionCreateAndPerform(mol, gMolActionDeleteImpropers, ig);
- IntGroupRelease(ig);
- return self;
+ 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:
- * assign_residue(group, res) -> Molecule
+ * find_conflicts(limit[, group1[, group2 [, ignore_exclusion]]]) -> [[n1, n2], [n3, n4], ...]
*
- * Assign the specified atoms as the given residue. res can either be an integer, "resname"
- * or "resname.resno". When the residue number is not specified, the residue number of
- * the first atom in the group is used.
- * This operation is undoable.
+ * 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_AssignResidue(VALUE self, VALUE range, VALUE res)
+s_Molecule_FindConflicts(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- IntGroup *ig;
- char *p, *pp, buf[16];
- Int resid, n;
- Atom *ap;
+ 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);
- /* Parse the argument res */
- if (FIXNUM_P(res)) {
- /* We can assume Fixnum here because Bignum is non-realistic as residue numbers */
- resid = NUM2INT(res);
- buf[0] = 0;
- } else {
- p = StringValuePtr(res);
- pp = strchr(p, '.');
- if (pp != NULL) {
- resid = atoi(pp + 1);
- n = pp - p;
- } else {
- resid = -1;
- n = strlen(p);
+ 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);
}
- if (n > sizeof buf - 1)
- n = sizeof buf - 1;
- strncpy(buf, p, n);
- buf[n] = 0;
- }
- ig = s_Molecule_AtomGroupFromValue(self, range);
- if (ig == NULL || IntGroupGetCount(ig) == 0)
- return Qnil;
-
- if (resid < 0) {
- /* Use the residue number of the first specified atom */
- n = IntGroupGetNthPoint(ig, 0);
- if (n >= mol->natoms)
- rb_raise(rb_eMolbyError, "Atom index (%d) out of range", n);
- ap = ATOM_AT_INDEX(mol->atoms, n);
- resid = ap->resSeq;
+ exinfo = mol->arena->exinfo; /* May be NULL */
+ exlist = mol->arena->exlist;
+ } else {
+ exinfo = NULL;
+ exlist = NULL;
}
- /* Change the residue number */
- MolActionCreateAndPerform(mol, gMolActionChangeResidueNumber, ig, resid);
- /* Change the residue name if necessary */
- if (buf[0] != 0) {
- /* Int seqs[2];
- seqs[0] = resid;
- seqs[1] = kInvalidIndex; */
- MolActionCreateAndPerform(mol, gMolActionChangeResidueNames, 1, &resid, 4, buf);
+ 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);
+ }
+ }
+ }
}
- IntGroupRelease(ig);
- return self;
-}
-
-/*
- * call-seq:
- * offset_residue(group, offset) -> Molecule
- *
- * Offset the residue number of the specified atoms. If any of the residue number gets
- * negative, then exception is thrown.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_OffsetResidue(VALUE self, VALUE range, VALUE offset)
-{
- Molecule *mol;
- IntGroup *ig;
- int ofs, result;
- Data_Get_Struct(self, Molecule, mol);
- ig = s_Molecule_AtomGroupFromValue(self, range);
- ofs = NUM2INT(offset);
- result = MolActionCreateAndPerform(mol, gMolActionOffsetResidueNumbers, ig, ofs, -1);
- if (result > 0)
- rb_raise(rb_eMolbyError, "residue number of atom %d becomes negative", result - 1);
- IntGroupRelease(ig);
- return self;
-}
-
-/*
- * call-seq:
- * renumber_atoms(array) -> IntGroup
- *
- * Change the order of atoms so that the atoms specified in the array argument appear
- * in this order from the top of the molecule. The atoms that are not included in array
- * are placed after these atoms, and these atoms are returned as an intGroup.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_RenumberAtoms(VALUE self, VALUE array)
-{
- Molecule *mol;
- Int *new2old;
- IntGroup *ig;
- int i, n;
- VALUE *valp, retval;
- Data_Get_Struct(self, Molecule, mol);
- if (TYPE(array) != T_ARRAY)
- array = rb_funcall(array, rb_intern("to_a"), 0);
- n = RARRAY_LEN(array);
- valp = RARRAY_PTR(array);
- new2old = ALLOC_N(Int, n + 1);
- for (i = 0; i < n; i++)
- new2old[i] = s_Molecule_AtomIndexFromValue(mol, valp[i]);
- new2old[i] = kInvalidIndex;
- i = MolActionCreateAndPerform(mol, gMolActionRenumberAtoms, i, new2old);
- if (i == 1)
- rb_raise(rb_eMolbyError, "Atom index out of range");
- else if (i == 2)
- rb_raise(rb_eMolbyError, "Duplicate entry");
- else if (i == 3)
- rb_raise(rb_eMolbyError, "Internal inconsistency during atom renumbering");
- retval = IntGroup_Alloc(rb_cIntGroup);
- Data_Get_Struct(retval, IntGroup, ig);
- if (mol->natoms > n)
- IntGroup_RaiseIfError(IntGroupAdd(ig, n, mol->natoms - n));
- xfree(new2old);
- return retval;
+ IntGroupIteratorRelease(&iter2);
+ IntGroupIteratorRelease(&iter1);
+ IntGroupRelease(ig2);
+ IntGroupRelease(ig1);
+ rval = rb_ary_new2(npairs);
+ if (pairs != NULL) {
+ for (i = 0; i < npairs; i++) {
+ rb_ary_push(rval, rb_ary_new3(2, INT2NUM(pairs[i * 2]), INT2NUM(pairs[i * 2 + 1])));
+ }
+ free(pairs);
+ }
+ return rval;
}
/*
if (nbonds > 0) {
for (n1 = 0; n1 < nbonds; n1++)
rb_ary_push(aval, INT2NUM(bonds[n1 * 2 + 1]));
- free(bonds);
- }
- return aval;
-}
-
-/*
- * call-seq:
- * guess_bonds(limit = 1.2) -> Integer
- *
- * Create bonds between atoms that are within the threshold distance.
- * If limit is a positive number, the threshold distance is the sum of the vdw radii times limit.
- * If limit is a negative number, its absolute value is used for the threshold distance in angstrom.
- * If limit is not given, a default value of 1.2 is used.
- * The number of the newly created bonds is returned.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_GuessBonds(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE limval;
- double limit;
- Int nbonds, *bonds;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &limval);
- if (limval == Qnil)
- limit = 1.2;
- else
- limit = NUM2DBL(rb_Float(limval));
- MoleculeGuessBonds(mol, limit, &nbonds, &bonds);
- if (nbonds > 0) {
- MolActionCreateAndPerform(mol, gMolActionAddBonds, nbonds * 2, bonds, NULL);
- free(bonds);
- }
- return INT2NUM(nbonds);
-}
-
-/*
- * call-seq:
- * register_undo(script, *args)
- *
- * Register an undo operation with the current molecule.
- */
-static VALUE
-s_Molecule_RegisterUndo(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE script, args;
- MolAction *act;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "1*", &script, &args);
- act = MolActionNew(SCRIPT_ACTION("R"), StringValuePtr(script), args);
- MolActionCallback_registerUndo(mol, act);
- return script;
-}
-
-/*
- * call-seq:
- * undo_enabled? -> bool
- *
- * Returns true if undo is enabled for this molecule; otherwise no.
- */
-static VALUE
-s_Molecule_UndoEnabled(VALUE self)
-{
- Molecule *mol;
- Data_Get_Struct(self, Molecule, mol);
- if (MolActionCallback_isUndoRegistrationEnabled(mol))
- return Qtrue;
- else return Qfalse;
-}
-
-/*
- * call-seq:
- * undo_enabled = bool
- *
- * Enable or disable undo.
- */
-static VALUE
-s_Molecule_SetUndoEnabled(VALUE self, VALUE val)
-{
- Molecule *mol;
- Data_Get_Struct(self, Molecule, mol);
- MolActionCallback_setUndoRegistrationEnabled(mol, (val != Qfalse && val != Qnil));
- return val;
-}
-
-/*
- * call-seq:
- * selection -> IntGroup
- *
- * Returns the current selection.
- */
-static VALUE
-s_Molecule_Selection(VALUE self)
-{
- Molecule *mol;
- IntGroup *ig;
- VALUE val;
- Data_Get_Struct(self, Molecule, mol);
- if (mol != NULL && (ig = MoleculeGetSelection(mol)) != NULL) {
- ig = IntGroupNewFromIntGroup(ig); /* Duplicate, so that the change from GUI does not affect the value */
- val = ValueFromIntGroup(ig);
- IntGroupRelease(ig);
- } else {
- val = IntGroup_Alloc(rb_cIntGroup);
- }
- return val;
-}
-
-static VALUE
-s_Molecule_SetSelectionSub(VALUE self, VALUE val, int undoable)
-{
- Molecule *mol;
- IntGroup *ig;
- Data_Get_Struct(self, Molecule, mol);
- if (val == Qnil)
- ig = NULL;
- else
- ig = s_Molecule_AtomGroupFromValue(self, val);
- if (undoable)
- MolActionCreateAndPerform(mol, gMolActionSetSelection, ig);
- else
- MoleculeSetSelection(mol, ig);
- if (ig != NULL)
- IntGroupRelease(ig);
- return val;
-}
-
-/*
- * call-seq:
- * selection = IntGroup
- *
- * Set the current selection. The right-hand operand may be nil.
- * This operation is _not_ undoable. If you need undo, use set_undoable_selection instead.
- */
-static VALUE
-s_Molecule_SetSelection(VALUE self, VALUE val)
-{
- return s_Molecule_SetSelectionSub(self, val, 0);
-}
-
-/*
- * call-seq:
- * set_undoable_selection(IntGroup)
- *
- * Set the current selection with undo registration. The right-hand operand may be nil.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_SetUndoableSelection(VALUE self, VALUE val)
-{
- return s_Molecule_SetSelectionSub(self, val, 1);
-}
-
-/*
- * call-seq:
- * hidden_atoms -> IntGroup
- *
- * Returns the currently hidden atoms.
- */
-static VALUE
-s_Molecule_HiddenAtoms(VALUE self)
-{
- rb_raise(rb_eMolbyError, "set_hidden_atoms is now obsolete. Try using Molecule#is_atom_visible or AtomRef#hidden.");
- return Qnil; /* Not reached */
-/* Molecule *mol;
- IntGroup *ig;
- VALUE val;
- Data_Get_Struct(self, Molecule, mol);
- if (mol != NULL) {
- Atom *ap;
- int i;
- ig = IntGroupNew();
- for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
- if (ap->exflags & kAtomHiddenFlag)
- IntGroupAdd(ig, i, 1);
- }
- val = ValueFromIntGroup(ig);
- IntGroupRelease(ig);
- rb_obj_freeze(val);
- return val;
- } else return Qnil; */
+ free(bonds);
+ }
+ return aval;
}
/*
* call-seq:
- * set_hidden_atoms(IntGroup)
- * self.hidden_atoms = IntGroup
+ * guess_bonds(limit = 1.2) -> Integer
*
- * Hide the specified atoms. This operation is _not_ undoable.
+ * Create bonds between atoms that are within the threshold distance.
+ * If limit is a positive number, the threshold distance is the sum of the vdw radii times limit.
+ * If limit is a negative number, its absolute value is used for the threshold distance in angstrom.
+ * If limit is not given, a default value of 1.2 is used.
+ * The number of the newly created bonds is returned.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_SetHiddenAtoms(VALUE self, VALUE val)
+s_Molecule_GuessBonds(int argc, VALUE *argv, VALUE self)
{
- rb_raise(rb_eMolbyError, "set_hidden_atoms is now obsolete. Try using Molecule#is_atom_visible or AtomRef#hidden.");
- return Qnil; /* Not reached */
-/*
- Molecule *mol;
+ Molecule *mol;
+ VALUE limval;
+ double limit;
+ Int nbonds, *bonds;
Data_Get_Struct(self, Molecule, mol);
- if (mol != NULL) {
- Atom *ap;
- int i;
- IntGroup *ig;
- if (val == Qnil)
- ig = NULL;
- else
- ig = s_Molecule_AtomGroupFromValue(self, val);
- for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
- if (ig != NULL && IntGroupLookup(ig, i, NULL)) {
- ap->exflags |= kAtomHiddenFlag;
- } else {
- ap->exflags &= kAtomHiddenFlag;
- }
- }
- if (ig != NULL)
- IntGroupRelease(ig);
- MoleculeCallback_notifyModification(mol, 0);
+ rb_scan_args(argc, argv, "01", &limval);
+ if (limval == Qnil)
+ limit = 1.2;
+ else
+ limit = NUM2DBL(rb_Float(limval));
+ MoleculeGuessBonds(mol, limit, &nbonds, &bonds);
+ if (nbonds > 0) {
+ MolActionCreateAndPerform(mol, gMolActionAddBonds, nbonds * 2, bonds, NULL);
+ free(bonds);
}
- return val; */
+ return INT2NUM(nbonds);
}
+#pragma mark ------ Cell and Symmetry ------
+
/*
* call-seq:
- * select_frame(index)
- * frame = index
+ * cell -> [a, b, c, alpha, beta, gamma [, sig_a, sig_b, sig_c, sig_alpha, sig_beta, sig_gamma]]
*
- * Select the specified frame. If successful, returns true, otherwise returns false.
+ * Returns the unit cell parameters. If cell is not set, returns nil.
*/
static VALUE
-s_Molecule_SelectFrame(VALUE self, VALUE val)
+s_Molecule_Cell(VALUE self)
{
Molecule *mol;
- int ival = NUM2INT(val);
+ int i;
+ VALUE val;
Data_Get_Struct(self, Molecule, mol);
- ival = MoleculeSelectFrame(mol, ival, 1);
- if (ival >= 0)
- return Qtrue;
- else return Qfalse;
+ if (mol->cell == NULL)
+ return Qnil;
+ val = rb_ary_new2(6);
+ for (i = 0; i < 6; i++)
+ rb_ary_push(val, rb_float_new(mol->cell->cell[i]));
+ if (mol->cell->has_sigma) {
+ for (i = 0; i < 6; i++) {
+ rb_ary_push(val, rb_float_new(mol->cell->cellsigma[i]));
+ }
+ }
+ return val;
}
/*
* call-seq:
- * frame -> Integer
+ * cell = [a, b, c, alpha, beta, gamma [, sig_a, sig_b, sig_c, sig_alpha, sig_beta, sig_gamma]]
+ * set_cell([a, b, c, alpha, beta, gamma[, sig_a, sig_b, sig_c, sig_alpha, sig_beta, sig_gamma]], convert_coord = nil)
*
- * Get the current frame.
+ * Set the unit cell parameters. If the cell value is nil, then clear the current cell.
+ If the given argument has 12 or more members, then the second half of the parameters represents the sigma values.
+ This operation is undoable.
+ Convert_coord is a flag to specify that the coordinates should be transformed so that the fractional coordinates remain the same.
*/
static VALUE
-s_Molecule_Frame(VALUE self)
+s_Molecule_SetCell(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
+ VALUE val, cval;
+ int i, convert_coord, n;
+ double d[12];
Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(mol->cframe);
+ rb_scan_args(argc, argv, "11", &val, &cval);
+ if (val == Qnil) {
+ n = 0;
+ } else {
+ int len;
+ val = rb_ary_to_ary(val);
+ len = RARRAY_LEN(val);
+ if (len >= 12) {
+ n = 12;
+ } else if (len >= 6) {
+ n = 6;
+ } else rb_raise(rb_eMolbyError, "too few members for cell parameters (6 or 12 required)");
+ for (i = 0; i < n; i++)
+ d[i] = NUM2DBL(rb_Float((RARRAY_PTR(val))[i]));
+ }
+ convert_coord = (RTEST(cval) ? 1 : 0);
+ MolActionCreateAndPerform(mol, gMolActionSetCell, n, d, convert_coord);
+ return val;
}
/*
* call-seq:
- * nframes -> Integer
+ * box -> [avec, bvec, cvec, origin, flags]
*
- * Get the number of frames.
+ * Get the unit cell information in the form of a periodic bounding box.
+ * Avec, bvec, cvec, origin are Vector3D objects, and flags is a 3-member array of
+ * Integers which define whether the system is periodic along the axis.
+ * If no unit cell is defined, nil is returned.
*/
static VALUE
-s_Molecule_Nframes(VALUE self)
+s_Molecule_Box(VALUE self)
{
Molecule *mol;
+ VALUE v[5], val;
Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(MoleculeGetNumberOfFrames(mol));
+ if (mol == NULL || mol->cell == NULL)
+ return Qnil;
+ v[0] = ValueFromVector(&(mol->cell->axes[0]));
+ v[1] = ValueFromVector(&(mol->cell->axes[1]));
+ v[2] = ValueFromVector(&(mol->cell->axes[2]));
+ v[3] = ValueFromVector(&(mol->cell->origin));
+ v[4] = rb_ary_new3(3, INT2NUM(mol->cell->flags[0]), INT2NUM(mol->cell->flags[1]), INT2NUM(mol->cell->flags[2]));
+ val = rb_ary_new4(5, v);
+ return val;
}
/*
* call-seq:
- * insert_frame(integer, coordinates = nil, cell_axes = nil) -> bool
- * insert_frames(intGroup = nil, coordinates = nil, cell_axes = nil) -> bool
+ * set_box(avec, bvec, cvec, origin = [0, 0, 0], flags = [1, 1, 1], convert_coordinates = nil)
+ * set_box(d, origin = [0, 0, 0])
+ * set_box
*
- * Insert new frames at the indices specified by the intGroup. If the first argument is
- * an integer, a single new frame is inserted at that index. If the first argument is
- * nil, a new frame is inserted at the last. If non-nil coordinates is given, it
- * should be an array of arrays of Vector3Ds, then those coordinates are set
- * to the new frame. Otherwise, the coordinates of current molecule are copied
- * to the new frame.
- * Returns an intGroup representing the inserted frames if successful, nil if not.
+ * Set the unit cell parameters. Avec, bvec, and cvec can be either a Vector3D or a number.
+ If it is a number, the x/y/z axis vector is multiplied with the given number and used
+ as the box vector.
+ Flags, if present, is a 3-member array of Integers defining whether the system is
+ periodic along the axis.
+ If convert_coordinates is true, then the coordinates are converted so that the fractional coordinates remain the same.
+ In the second form, an isotropic box with cell-length d is set.
+ In the third form, the existing box is cleared.
+ Note: the sigma of the cell parameters is not cleared unless the periodic box itself is cleared.
*/
static VALUE
-s_Molecule_InsertFrames(int argc, VALUE *argv, VALUE self)
+s_Molecule_SetBox(VALUE self, VALUE aval)
{
- VALUE val, coords, cells;
Molecule *mol;
- IntGroup *ig;
- int count, ival, i, j, len, len_c, len2, nframes;
- VALUE *ptr, *ptr2;
- Vector *vp, *vp2;
+ VALUE v[6];
+ static Vector ax[3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+ Vector vv[3];
+ Vector origin = {0, 0, 0};
+ char flags[3];
+ Double d;
+ int i, convertCoordinates = 0;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "12", &val, &coords, &cells);
- if (coords != Qnil) {
- if (TYPE(coords) != T_ARRAY)
- rb_raise(rb_eTypeError, "the coordinates should be given as an array of Vector3D");
- len = RARRAY_LEN(coords);
- } else len = 0;
- if (cells != Qnil) {
- if (mol->cell == NULL)
- rb_raise(rb_eTypeError, "the unit cell is not defined but the cell axes are given");
- if (TYPE(cells) != T_ARRAY)
- rb_raise(rb_eTypeError, "the cell axes should be given as an array of Vector3D");
- len_c = RARRAY_LEN(cells);
- } else len_c = 0;
- count = (len > len_c ? len : len_c); /* May be zero; will be updated later */
- nframes = MoleculeGetNumberOfFrames(mol);
- if (val == Qnil) {
- ig = IntGroupNewWithPoints(nframes, (count > 0 ? count : 1), -1);
- val = ValueFromIntGroup(ig);
- } else {
- ig = IntGroupFromValue(val);
+ if (aval == Qnil) {
+ MolActionCreateAndPerform(mol, gMolActionClearBox);
+ return self;
}
- count = IntGroupGetCount(ig); /* Count is updated here */
- vp = ALLOC_N(Vector, mol->natoms * count);
- if (cells != Qnil)
- vp2 = ALLOC_N(Vector, 4 * count);
- else vp2 = NULL;
- if (len > 0) {
- if (len < count)
- rb_raise(rb_eMolbyError, "the coordinates should contain no less than %d arrays (for frames)", count);
- ptr = RARRAY_PTR(coords);
- for (i = 0; i < count; i++) {
- if (TYPE(ptr[i]) != T_ARRAY)
- rb_raise(rb_eTypeError, "the coordinate array contains non-array object at index %d", i);
- len2 = RARRAY_LEN(ptr[i]);
- if (len2 < mol->natoms)
- rb_raise(rb_eMolbyError, "the array at index %d contains less than %d elements", i, mol->natoms);
- ptr2 = RARRAY_PTR(ptr[i]);
- for (j = 0; j < mol->natoms; j++)
- VectorFromValue(ptr2[j], &vp[i * mol->natoms + j]);
- }
- } else {
- Atom *ap;
- for (i = 0; i < count; i++) {
- for (j = 0, ap = mol->atoms; j < mol->natoms; j++, ap = ATOM_NEXT(ap)) {
- vp[i * mol->natoms + j] = ap->r;
- }
- }
+ aval = rb_ary_to_ary(aval);
+ for (i = 0; i < 6; i++) {
+ if (i < RARRAY_LEN(aval))
+ v[i] = (RARRAY_PTR(aval))[i];
+ else v[i] = Qnil;
}
- if (len_c > 0) {
- if (len_c < count)
- rb_raise(rb_eMolbyError, "the cell vectors should contain no less than %d arrays (for frames)", count);
- ptr = RARRAY_PTR(cells);
- for (i = 0; i < count; i++) {
- if (TYPE(ptr[i]) != T_ARRAY)
- rb_raise(rb_eTypeError, "the cell parameter array contains non-array object at index %d", i);
- len2 = RARRAY_LEN(ptr[i]);
- if (len2 < 4)
- rb_raise(rb_eMolbyError, "the cell parameter should contain 4 vectors");
- ptr2 = RARRAY_PTR(ptr[i]);
- for (j = 0; j < 4; j++)
- VectorFromValue(ptr2[j], &vp2[i * 4 + j]);
- }
+ if (v[0] == Qnil) {
+ MolActionCreateAndPerform(mol, gMolActionClearBox);
+ return self;
}
- ival = MolActionCreateAndPerform(mol, gMolActionInsertFrames, ig, mol->natoms * count, vp, (vp2 != NULL ? 4 * count : 0), vp2);
- IntGroupRelease(ig);
- xfree(vp);
- if (vp2 != NULL)
- xfree(vp2);
- return (ival >= 0 ? val : Qnil);
-}
-
-/*
- * call-seq:
- * create_frame(coordinates = nil) -> Integer
- * create_frames(coordinates = nil) -> Integer
- *
- * Same as molecule.insert_frames(nil, coordinates).
- */
-static VALUE
-s_Molecule_CreateFrames(int argc, VALUE *argv, VALUE self)
-{
- VALUE vals[3];
- rb_scan_args(argc, argv, "02", &vals[1], &vals[2]);
- vals[0] = Qnil;
- return s_Molecule_InsertFrames(3, vals, self);
+ if ((v[1] == Qnil || v[2] == Qnil) && rb_obj_is_kind_of(v[0], rb_cNumeric)) {
+ d = NUM2DBL(rb_Float(v[0]));
+ for (i = 0; i < 3; i++)
+ VecScale(vv[i], ax[i], d);
+ if (v[1] != Qnil)
+ VectorFromValue(v[1], &origin);
+ flags[0] = flags[1] = flags[2] = 1;
+ } else {
+ for (i = 0; i < 3; i++) {
+ if (v[i] == Qnil) {
+ VecZero(vv[i]);
+ } else if (rb_obj_is_kind_of(v[i], rb_cNumeric)) {
+ d = NUM2DBL(rb_Float(v[i]));
+ VecScale(vv[i], ax[i], d);
+ } else {
+ VectorFromValue(v[i], &vv[i]);
+ }
+ flags[i] = (VecLength2(vv[i]) > 0.0);
+ }
+ if (v[3] != Qnil)
+ VectorFromValue(v[3], &origin);
+ if (v[4] != Qnil) {
+ for (i = 0; i < 3; i++) {
+ VALUE val = Ruby_ObjectAtIndex(v[4], i);
+ flags[i] = (NUM2INT(rb_Integer(val)) != 0);
+ }
+ }
+ if (RTEST(v[5]))
+ convertCoordinates = 1;
+ }
+ MolActionCreateAndPerform(mol, gMolActionSetBox, &(vv[0]), &(vv[1]), &(vv[2]), &origin, (flags[0] * 4 + flags[1] * 2 + flags[2]), convertCoordinates);
+ return self;
}
/*
* call-seq:
- * remove_frames(IntGroup, wantCoordinates = false)
+ * cell_periodicity -> [n1, n2, n3]
*
- * Remove the frames at group. If wantsCoordinates is false (default), returns true if successful
- * and nil otherwise. If wantsCoordinates is true, an array of arrays of the coordinates in the
- * removed frames is returned if operation is successful.
+ * Get flags denoting whether the cell is periodic along the a/b/c axes. If the cell is not defined
+ * nil is returned.
*/
static VALUE
-s_Molecule_RemoveFrames(int argc, VALUE *argv, VALUE self)
+s_Molecule_CellPeriodicity(VALUE self)
{
- VALUE val, flag;
- VALUE retval;
Molecule *mol;
- IntGroup *ig;
- int count;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "11", &val, &flag);
- ig = IntGroupFromValue(val);
- count = IntGroupGetCount(ig);
- if (RTEST(flag)) {
- /* Create return value before removing frames */
- VALUE coords;
- int i, j, n;
- Atom *ap;
- Vector v;
- retval = rb_ary_new2(count);
- for (i = 0; i < count; i++) {
- n = IntGroupGetNthPoint(ig, i);
- coords = rb_ary_new2(mol->natoms);
- for (j = 0, ap = mol->atoms; j < mol->natoms; j++, ap = ATOM_NEXT(ap)) {
- if (n < ap->nframes && n != mol->cframe)
- v = ap->frames[n];
- else v = ap->r;
- rb_ary_push(coords, ValueFromVector(&v));
- }
- rb_ary_push(retval, coords);
- }
- } else retval = Qtrue;
- if (MolActionCreateAndPerform(mol, gMolActionRemoveFrames, ig) >= 0)
- return retval;
- else return Qnil;
+ if (mol->cell == NULL)
+ return Qnil;
+ return rb_ary_new3(3, INT2FIX((int)mol->cell->flags[0]), INT2FIX((int)mol->cell->flags[1]), INT2FIX((int)mol->cell->flags[2]));
}
/*
* call-seq:
- * each_frame {|n| ...}
+ * self.cell_periodicity = [n1, n2, n3] or Integer or nil
+ * set_cell_periodicity = [n1, n2, n3] or Integer or nil
*
- * Set the frame number from 0 to nframes-1 and execute the block. The block argument is
- * the frame number. After completion, the original frame number is restored.
+ * Set whether the cell is periodic along the a/b/c axes. If an integer is given as an argument,
+ * its bits 2/1/0 (from the lowest) correspond to the a/b/c axes. Nil is equivalent to [0, 0, 0].
+ * If cell is not defined, exception is raised.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_EachFrame(VALUE self)
+s_Molecule_SetCellPeriodicity(VALUE self, VALUE arg)
{
- int i, cframe, nframes;
Molecule *mol;
+ Int flag;
Data_Get_Struct(self, Molecule, mol);
- cframe = mol->cframe;
- nframes = MoleculeGetNumberOfFrames(mol);
- if (nframes > 0) {
- for (i = 0; i < nframes; i++) {
- MoleculeSelectFrame(mol, i, 1);
- rb_yield(INT2NUM(i));
+ if (mol->cell == NULL)
+ rb_raise(rb_eMolbyError, "periodic cell is not defined");
+ if (arg == Qnil)
+ flag = 0;
+ else if (rb_obj_is_kind_of(arg, rb_cNumeric))
+ flag = NUM2INT(rb_Integer(arg));
+ else {
+ Int i;
+ VALUE arg0;
+ arg = rb_ary_to_ary(arg);
+ flag = 0;
+ for (i = 0; i < 3 && i < RARRAY_LEN(arg); i++) {
+ arg0 = RARRAY_PTR(arg)[i];
+ if (arg0 != Qnil && arg0 != Qfalse && arg0 != INT2FIX(0))
+ flag |= (1 << (2 - i));
}
- MoleculeSelectFrame(mol, cframe, 1);
}
- return self;
+ MolActionCreateAndPerform(mol, gMolActionSetCellPeriodicity, flag);
+ return arg;
}
/*
* call-seq:
- * get_coord_from_frame(index, group = nil)
+ * cell_flexibility -> bool
*
- * Copy the coordinates from the indicated frame. If group is specified, only the specified atoms
- * are modified. Third argument (cflag) is now obsolete (it used to specify whether the cell parameters are to be
- * copied; now they are always copied)
+ * Returns the unit cell is flexible or not
*/
static VALUE
-s_Molecule_GetCoordFromFrame(int argc, VALUE *argv, VALUE self)
+s_Molecule_CellFlexibility(VALUE self)
{
- Molecule *mol;
- VALUE ival, gval, cval;
- Int index, i, j, n, nn;
- IntGroup *ig;
- IntGroupIterator iter;
- Atom *ap;
- Vector *vp;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "12", &ival, &gval, &cval);
- if (argc == 3)
- rb_warn("The 3rd argument to get_coord_from_frame() is now obsolete");
- index = NUM2INT(rb_Integer(ival));
- if (index < 0 || index >= (n = MoleculeGetNumberOfFrames(mol))) {
- if (n == 0)
- rb_raise(rb_eMolbyError, "No frame is present");
- else
- rb_raise(rb_eMolbyError, "Frame index (%d) out of range (should be 0..%d)", index, n - 1);
- }
- if (gval == Qnil) {
- ig = IntGroupNewWithPoints(0, mol->natoms, -1);
- } else {
- ig = s_Molecule_AtomGroupFromValue(self, gval);
- }
- n = IntGroupGetCount(ig);
- if (n > 0) {
- vp = (Vector *)calloc(sizeof(Vector), n);
- IntGroupIteratorInit(ig, &iter);
- j = 0;
- nn = 0;
- while ((i = IntGroupIteratorNext(&iter)) >= 0) {
- ap = ATOM_AT_INDEX(mol->atoms, i);
- if (index < ap->nframes) {
- vp[j] = ap->frames[index];
- nn++;
- } else {
- vp[j] = ap->r;
- }
- j++;
- }
- if (nn > 0)
- MolActionCreateAndPerform(mol, gMolActionSetAtomPositions, ig, n, vp);
- free(vp);
- if (mol->cell != NULL && mol->frame_cells != NULL && index < mol->nframe_cells) {
- vp = mol->frame_cells + index * 4;
- MolActionCreateAndPerform(mol, gMolActionSetBox, vp, vp + 1, vp + 2, vp + 3, -1, 0);
- }
- IntGroupIteratorRelease(&iter);
- }
- /* Copy the extra properties */
- IntGroupRelease(ig);
- for (i = 0; i < mol->nmolprops; i++) {
- Double *dp = (Double *)malloc(sizeof(Double));
- ig = IntGroupNew();
- IntGroupAdd(ig, mol->cframe, 1);
- *dp = mol->molprops[i].propvals[index];
- MolActionCreateAndPerform(mol, gMolActionSetProperty, i, ig, 1, dp);
- free(dp);
- IntGroupRelease(ig);
- }
-
- return self;
+ rb_warn("cell_flexibility is obsolete (unit cell is always frame dependent)");
+ return Qtrue;
+ /* Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->cell == NULL)
+ return Qfalse;
+ if (mol->useFlexibleCell)
+ return Qtrue;
+ else return Qfalse; */
}
/*
* call-seq:
- * reorder_frames(old_indices)
+ * self.cell_flexibility = bool
+ * set_cell_flexibility(bool)
*
- * Reorder the frames. The argument is an array of integers that specify the 'old'
- * frame numbers. Thus, if the argument is [2,0,1], then the new frames 0/1/2 are the
- * same as the old frames 2/0/1, respectively.
- * The argument must have the same number of integers as the number of frames.
+ * Change the unit cell is flexible or not
*/
static VALUE
-s_Molecule_ReorderFrames(VALUE self, VALUE aval)
+s_Molecule_SetCellFlexibility(VALUE self, VALUE arg)
{
- Molecule *mol;
- Int *ip, *ip2, i, n, nframes;
- Data_Get_Struct(self, Molecule, mol);
- aval = rb_ary_to_ary(aval);
- nframes = MoleculeGetNumberOfFrames(mol);
- if (RARRAY_LEN(aval) != nframes)
- rb_raise(rb_eMolbyError, "The argument must have the same number of integers as the number of frames");
- ip2 = (Int *)calloc(sizeof(Int), nframes);
- ip = (Int *)calloc(sizeof(Int), nframes);
- for (i = 0; i < nframes; i++) {
- n = NUM2INT(rb_Integer(RARRAY_PTR(aval)[i]));
- if (n < 0 || n >= nframes || ip2[n] != 0) {
- free(ip2);
- free(ip);
- if (n < 0 || n >= nframes)
- rb_raise(rb_eMolbyError, "The argument (%d) is out of range", n);
- else
- rb_raise(rb_eMolbyError, "The argument has duplicated entry (%d)", n);
- }
- ip2[n] = 1;
- ip[i] = n;
- }
- free(ip2);
- MolActionCreateAndPerform(mol, gMolActionReorderFrames, nframes, ip);
- free(ip);
+ rb_warn("set_cell_flexibility is obsolete (unit cell is always frame dependent)");
return self;
+ /* Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ MolActionCreateAndPerform(mol, gMolActionSetCellFlexibility, RTEST(arg) != 0);
+ return self; */
}
-
+
/*
* call-seq:
- * set_atom_attr(index, key, value)
+ * cell_transform -> Transform
*
- * Set the atom attribute for the specified atom.
- * This operation is undoable.
+ * Get the transform matrix that converts internal coordinates to cartesian coordinates.
+ * If cell is not defined, nil is returned.
*/
static VALUE
-s_Molecule_SetAtomAttr(VALUE self, VALUE idx, VALUE key, VALUE val)
+s_Molecule_CellTransform(VALUE self)
{
- Molecule *mol;
- VALUE aref, oldval;
+ Molecule *mol;
Data_Get_Struct(self, Molecule, mol);
- aref = ValueFromMoleculeAndIndex(mol, s_Molecule_AtomIndexFromValue(mol, idx));
- oldval = s_AtomRef_GetAttr(aref, key);
- if (val == Qundef)
- return oldval;
- s_AtomRef_SetAttr(aref, key, val);
- return val;
+ if (mol == NULL || mol->cell == NULL)
+ return Qnil;
+ return ValueFromTransform(&(mol->cell->tr));
}
/*
* call-seq:
- * get_atom_attr(index, key)
+ * symmetry -> Array of Transforms
+ * symmetries -> Array of Transforms
*
- * Get the atom attribute for the specified atom.
+ * Get the currently defined symmetry operations. If no symmetry operation is defined,
+ * returns an empty array.
*/
static VALUE
-s_Molecule_GetAtomAttr(VALUE self, VALUE idx, VALUE key)
+s_Molecule_Symmetry(VALUE self)
{
- return s_Molecule_SetAtomAttr(self, idx, key, Qundef);
+ Molecule *mol;
+ VALUE val;
+ int i;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->nsyms <= 0)
+ return rb_ary_new();
+ val = rb_ary_new2(mol->nsyms);
+ for (i = 0; i < mol->nsyms; i++) {
+ rb_ary_push(val, ValueFromTransform(&mol->syms[i]));
+ }
+ return val;
}
/*
* call-seq:
- * fragment(n1, *exatoms) -> IntGroup
- * fragment(group, *exatoms) -> IntGroup
+ * nsymmetries -> Integer
*
- * Get the fragment including the atom n1 or the atom group. If additional arguments are given,
- * those atoms will not be counted during the search.
+ * Get the number of currently defined symmetry operations.
*/
static VALUE
-s_Molecule_Fragment(int argc, VALUE *argv, VALUE self)
+s_Molecule_Nsymmetries(VALUE self)
{
Molecule *mol;
- IntGroup *baseg, *ig, *exatoms;
- int n;
- volatile VALUE nval, exval;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "1*", &nval, &exval);
- if (rb_obj_is_kind_of(nval, rb_cNumeric) || rb_obj_is_kind_of(nval, rb_cString)) {
- baseg = NULL;
- n = NUM2INT(s_Molecule_AtomIndex(self, nval));
- } else {
- baseg = s_Molecule_AtomGroupFromValue(self, nval);
- }
- if (RARRAY_LEN(exval) == 0) {
- exatoms = NULL;
- } else {
- exval = s_Molecule_AtomGroup(RARRAY_LEN(exval), RARRAY_PTR(exval), self);
- Data_Get_Struct(exval, IntGroup, exatoms);
- }
- if (baseg == NULL) {
- ig = MoleculeFragmentExcludingAtomGroup(mol, n, exatoms);
- } else {
- IntGroupIterator iter;
- IntGroupIteratorInit(baseg, &iter);
- if ((n = IntGroupIteratorNext(&iter)) < 0) {
- ig = IntGroupNew();
- } else {
- ig = MoleculeFragmentExcludingAtomGroup(mol, n, exatoms);
- if (ig != NULL) {
- while ((n = IntGroupIteratorNext(&iter)) >= 0) {
- IntGroup *subg;
- subg = MoleculeFragmentExcludingAtomGroup(mol, n, exatoms);
- if (subg != NULL) {
- IntGroupAddIntGroup(ig, subg);
- IntGroupRelease(subg);
- }
- }
- }
- }
- IntGroupIteratorRelease(&iter);
- IntGroupRelease(baseg);
- }
- if (ig == NULL)
- rb_raise(rb_eMolbyError, "invalid specification of molecular fragment");
- nval = ValueFromIntGroup(ig);
- IntGroupRelease(ig);
- return nval;
+ return INT2NUM(mol->nsyms);
}
/*
* call-seq:
- * fragments(exclude = nil)
+ * add_symmetry(Transform) -> Integer
*
- * Returns the fragments as an array of IntGroups.
- * If exclude is given (as an array or an IntGroup), then those atoms are excluded
- * in defining the fragment.
+ * Add a new symmetry operation. If no symmetry operation is defined and the
+ * given argument is not an identity transform, then also add an identity
+ * transform at the index 0.
+ * Returns the total number of symmetries after operation.
*/
static VALUE
-s_Molecule_Fragments(int argc, VALUE *argv, VALUE self)
+s_Molecule_AddSymmetry(VALUE self, VALUE trans)
{
Molecule *mol;
- IntGroup *ag, *fg, *eg;
- VALUE gval, exval, retval;
+ Transform tr;
Data_Get_Struct(self, Molecule, mol);
- if (mol == NULL)
- return Qnil;
- if (mol->natoms == 0)
- return rb_ary_new();
- rb_scan_args(argc, argv, "01", &exval);
- if (exval == Qnil)
- eg = NULL;
- else
- eg = IntGroupFromValue(exval);
- ag = IntGroupNewWithPoints(0, mol->natoms, -1);
- if (eg != NULL)
- IntGroupRemoveIntGroup(ag, eg);
- retval = rb_ary_new();
- while (IntGroupGetCount(ag) > 0) {
- int n = IntGroupGetNthPoint(ag, 0);
- fg = MoleculeFragmentExcludingAtomGroup(mol, n, eg);
- if (fg == NULL)
- rb_raise(rb_eMolbyError, "internal error during each_fragment");
- gval = ValueFromIntGroup(fg);
- rb_ary_push(retval, gval);
- IntGroupRemoveIntGroup(ag, fg);
- IntGroupRelease(fg);
- }
- IntGroupRelease(ag);
- if (eg != NULL)
- IntGroupRelease(eg);
- return retval;
+ TransformFromValue(trans, &tr);
+ MolActionCreateAndPerform(mol, gMolActionAddSymmetryOperation, &tr);
+ return INT2NUM(mol->nsyms);
}
/*
* call-seq:
- * each_fragment(exclude = nil) {|group| ...}
+ * remove_symmetry(count = nil) -> Integer
+ * remove_symmetries(count = nil) -> Integer
*
- * Execute the block, with the IntGroup object for each fragment as the argument.
- * Atoms or bonds should not be added or removed during the execution of the block.
- * If exclude is given (as an array or an IntGroup), then those atoms are excluded
- * in defining the fragment.
+ * Remove the specified number of symmetry operations. The last added ones are removed
+ * first. If count is nil, then all symmetry operations are removed. Returns the
+ * number of leftover symmetries.
*/
static VALUE
-s_Molecule_EachFragment(int argc, VALUE *argv, VALUE self)
+s_Molecule_RemoveSymmetry(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- IntGroup *ag, *fg, *eg;
- VALUE gval, exval;
+ VALUE cval;
+ int i, n;
Data_Get_Struct(self, Molecule, mol);
- if (mol == NULL || mol->natoms == 0)
- return self;
- rb_scan_args(argc, argv, "01", &exval);
- if (exval == Qnil)
- eg = NULL;
- else
- eg = IntGroupFromValue(exval);
- ag = IntGroupNewWithPoints(0, mol->natoms, -1);
- if (eg != NULL)
- IntGroupRemoveIntGroup(ag, eg);
- while (IntGroupGetCount(ag) > 0) {
- int n = IntGroupGetNthPoint(ag, 0);
- fg = MoleculeFragmentExcludingAtomGroup(mol, n, eg);
- if (fg == NULL)
- rb_raise(rb_eMolbyError, "internal error during each_fragment");
- gval = ValueFromIntGroup(fg);
- rb_yield(gval);
- IntGroupRemoveIntGroup(ag, fg);
- IntGroupRelease(fg);
+ rb_scan_args(argc, argv, "01", &cval);
+ if (cval == Qnil)
+ n = mol->nsyms - 1;
+ else {
+ n = NUM2INT(rb_Integer(cval));
+ if (n < 0 || n > mol->nsyms)
+ rb_raise(rb_eMolbyError, "the given count of symops is out of range");
+ if (n == mol->nsyms)
+ n = mol->nsyms - 1;
}
- IntGroupRelease(ag);
- if (eg != NULL)
- IntGroupRelease(eg);
- return self;
+ for (i = 0; i < n; i++)
+ MolActionCreateAndPerform(mol, gMolActionDeleteSymmetryOperation);
+ return INT2NUM(mol->nsyms);
}
/*
* call-seq:
- * detachable?(group) -> [n1, n2]
+ * wrap_unit_cell(group) -> Vector3D
*
- * Check whether the group is 'detachable', i.e. the group is bound to the rest
- * of the molecule via only one bond. If it is, then the indices of the atoms
- * belonging to the bond is returned, the first element being the atom included
- * in the fragment. Otherwise, Qnil is returned.
+ * 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_Detachable_P(VALUE self, VALUE gval)
+s_Molecule_WrapUnitCell(VALUE self, VALUE gval)
{
- Molecule *mol;
+ Molecule *mol;
IntGroup *ig;
- int n1, n2;
- VALUE retval;
+ 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);
- if (MoleculeIsFragmentDetachable(mol, ig, &n1, &n2)) {
- retval = rb_ary_new3(2, INT2NUM(n1), INT2NUM(n2));
- } else retval = Qnil;
+ 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 retval;
+ return ValueFromVector(&dv);
}
/*
* call-seq:
- * bonds_on_border(group = selection) -> Array of Array of two Integers
+ * expand_by_symmetry(group, sym, dx=0, dy=0, dz=0, allow_overlap = false) -> Array
*
- * 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.
+ * Expand the specified part of the molecule by the given symmetry operation.
+ * Returns the array of atom indices corresponding to the expanded atoms.
+ * If allow_overlap is true, then new atoms are created even when the
+ * coordinates coincide with the some other atom (special position) of the
+ * same element; otherwise, such atom will not be created and the index of the
+ * existing atom is given in the returned array.
*/
static VALUE
-s_Molecule_BondsOnBorder(int argc, VALUE *argv, VALUE self)
+s_Molecule_ExpandBySymmetry(int argc, VALUE *argv, VALUE self)
{
- Molecule *mol;
- IntGroup *ig, *bg;
- VALUE gval, retval;
+ Molecule *mol;
+ VALUE gval, sval, xval, yval, zval, rval, oval;
+ IntGroup *ig;
+ Int n[4], allow_overlap;
+ Int natoms;
+ Int nidx, *idx;
+
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);
+ rb_scan_args(argc, argv, "24", &gval, &sval, &xval, &yval, &zval, &oval);
+ 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)));
+ allow_overlap = (RTEST(oval) ? 1 : 0);
+ 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], allow_overlap, &nidx, &idx);
+
+ rval = rb_ary_new2(nidx);
+ while (--nidx >= 0) {
+ rb_ary_store(rval, nidx, INT2NUM(idx[nidx]));
}
- IntGroupRelease(bg);
- IntGroupRelease(ig);
- return retval;
+ /* 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;
}
+#pragma mark ------ Transforms ------
+
/*
* call-seq:
* translate(vec, group = nil) -> Molecule
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);
+ // MoleculeTranslate(mol, &v, ig);
MolActionCreateAndPerform(mol, gMolActionTranslateAtoms, &v, ig);
if (ig != NULL)
IntGroupRelease(ig);
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); */
+ /* 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)
+/*
+ * 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)
{
- 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;
+ 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:
- * center_of_mass(group = nil) -> Vector3D
+ * symop_for_transform(transform, is_cartesian = nil) -> [sym, dx, dy, dz]
*
- * Calculate the center of mass for the given set of atoms. The argument
- * group is null, then all atoms are considered.
+ * 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_CenterOfMass(int argc, VALUE *argv, VALUE self)
+s_Molecule_SymopForTransform(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE group;
- IntGroup *ig;
- Vector v;
+ VALUE tval, fval;
+ Symop symop;
+ Transform tr;
+ int n;
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);
+ 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 */
+ }
}
+#pragma mark ------ Frames ------
+
/*
* call-seq:
- * centralize(group = nil) -> self
+ * select_frame(index)
+ * frame = index
*
- * 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).
+ * Select the specified frame. If successful, returns true, otherwise returns false.
*/
static VALUE
-s_Molecule_Centralize(int argc, VALUE *argv, VALUE self)
+s_Molecule_SelectFrame(VALUE self, VALUE val)
{
Molecule *mol;
- VALUE group;
- IntGroup *ig;
- Vector v;
+ int ival = NUM2INT(val);
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;
+ ival = MoleculeSelectFrame(mol, ival, 1);
+ if (ival >= 0)
+ return Qtrue;
+ else return Qfalse;
}
/*
* call-seq:
- * bounds(group = nil) -> [min, max]
+ * frame -> Integer
*
- * Calculate the boundary. The return value is an array of two Vector3D objects.
+ * Get the current frame.
*/
static VALUE
-s_Molecule_Bounds(int argc, VALUE *argv, VALUE self)
+s_Molecule_Frame(VALUE self)
{
Molecule *mol;
- VALUE group;
+ Data_Get_Struct(self, Molecule, mol);
+ return INT2NUM(mol->cframe);
+}
+
+/*
+ * call-seq:
+ * nframes -> Integer
+ *
+ * Get the number of frames.
+ */
+static VALUE
+s_Molecule_Nframes(VALUE self)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ return INT2NUM(MoleculeGetNumberOfFrames(mol));
+}
+
+/*
+ * call-seq:
+ * insert_frame(integer, coordinates = nil, cell_axes = nil) -> bool
+ * insert_frames(intGroup = nil, coordinates = nil, cell_axes = nil) -> bool
+ *
+ * Insert new frames at the indices specified by the intGroup. If the first argument is
+ * an integer, a single new frame is inserted at that index. If the first argument is
+ * nil, a new frame is inserted at the last. If non-nil coordinates is given, it
+ * should be an array of arrays of Vector3Ds, then those coordinates are set
+ * to the new frame. Otherwise, the coordinates of current molecule are copied
+ * to the new frame.
+ * Returns an intGroup representing the inserted frames if successful, nil if not.
+ */
+static VALUE
+s_Molecule_InsertFrames(int argc, VALUE *argv, VALUE self)
+{
+ VALUE val, coords, cells;
+ Molecule *mol;
IntGroup *ig;
- Vector vmin, vmax;
- int n;
- Atom *ap;
+ int count, ival, i, j, len, len_c, len2, nframes;
+ VALUE *ptr, *ptr2;
+ Vector *vp, *vp2;
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;
+ rb_scan_args(argc, argv, "12", &val, &coords, &cells);
+ if (coords != Qnil) {
+ if (TYPE(coords) != T_ARRAY)
+ rb_raise(rb_eTypeError, "the coordinates should be given as an array of Vector3D");
+ len = RARRAY_LEN(coords);
+ } else len = 0;
+ if (cells != Qnil) {
+ if (mol->cell == NULL)
+ rb_raise(rb_eTypeError, "the unit cell is not defined but the cell axes are given");
+ if (TYPE(cells) != T_ARRAY)
+ rb_raise(rb_eTypeError, "the cell axes should be given as an array of Vector3D");
+ len_c = RARRAY_LEN(cells);
+ } else len_c = 0;
+ count = (len > len_c ? len : len_c); /* May be zero; will be updated later */
+ nframes = MoleculeGetNumberOfFrames(mol);
+ if (val == Qnil) {
+ ig = IntGroupNewWithPoints(nframes, (count > 0 ? count : 1), -1);
+ val = ValueFromIntGroup(ig);
+ } else {
+ ig = IntGroupFromValue(val);
+ }
+ count = IntGroupGetCount(ig); /* Count is updated here */
+ vp = ALLOC_N(Vector, mol->natoms * count);
+ if (cells != Qnil)
+ vp2 = ALLOC_N(Vector, 4 * count);
+ else vp2 = NULL;
+ if (len > 0) {
+ if (len < count)
+ rb_raise(rb_eMolbyError, "the coordinates should contain no less than %d arrays (for frames)", count);
+ ptr = RARRAY_PTR(coords);
+ for (i = 0; i < count; i++) {
+ if (TYPE(ptr[i]) != T_ARRAY)
+ rb_raise(rb_eTypeError, "the coordinate array contains non-array object at index %d", i);
+ len2 = RARRAY_LEN(ptr[i]);
+ if (len2 < mol->natoms)
+ rb_raise(rb_eMolbyError, "the array at index %d contains less than %d elements", i, mol->natoms);
+ ptr2 = RARRAY_PTR(ptr[i]);
+ for (j = 0; j < mol->natoms; j++)
+ VectorFromValue(ptr2[j], &vp[i * mol->natoms + j]);
+ }
+ } else {
+ Atom *ap;
+ for (i = 0; i < count; i++) {
+ for (j = 0, ap = mol->atoms; j < mol->natoms; j++, ap = ATOM_NEXT(ap)) {
+ vp[i * mol->natoms + j] = ap->r;
+ }
+ }
+ }
+ if (len_c > 0) {
+ if (len_c < count)
+ rb_raise(rb_eMolbyError, "the cell vectors should contain no less than %d arrays (for frames)", count);
+ ptr = RARRAY_PTR(cells);
+ for (i = 0; i < count; i++) {
+ if (TYPE(ptr[i]) != T_ARRAY)
+ rb_raise(rb_eTypeError, "the cell parameter array contains non-array object at index %d", i);
+ len2 = RARRAY_LEN(ptr[i]);
+ if (len2 < 4)
+ rb_raise(rb_eMolbyError, "the cell parameter should contain 4 vectors");
+ ptr2 = RARRAY_PTR(ptr[i]);
+ for (j = 0; j < 4; j++)
+ VectorFromValue(ptr2[j], &vp2[i * 4 + j]);
+ }
}
- return rb_ary_new3(2, ValueFromVector(&vmin), ValueFromVector(&vmax));
+ ival = MolActionCreateAndPerform(mol, gMolActionInsertFrames, ig, mol->natoms * count, vp, (vp2 != NULL ? 4 * count : 0), vp2);
+ IntGroupRelease(ig);
+ xfree(vp);
+ if (vp2 != NULL)
+ xfree(vp2);
+ return (ival >= 0 ? val : Qnil);
}
-/* Get atom position or a vector */
-static void
-s_Molecule_GetVectorFromArg(Molecule *mol, VALUE val, Vector *vp)
+/*
+ * call-seq:
+ * create_frame(coordinates = nil) -> Integer
+ * create_frames(coordinates = nil) -> Integer
+ *
+ * Same as molecule.insert_frames(nil, coordinates).
+ */
+static VALUE
+s_Molecule_CreateFrames(int argc, VALUE *argv, VALUE self)
{
- 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);
- }
+ VALUE vals[3];
+ rb_scan_args(argc, argv, "02", &vals[1], &vals[2]);
+ vals[0] = Qnil;
+ return s_Molecule_InsertFrames(3, vals, self);
}
/*
* call-seq:
- * measure_bond(n1, n2) -> Float
+ * remove_frames(IntGroup, wantCoordinates = false)
*
- * 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.
+ * Remove the frames at group. If wantsCoordinates is false (default), returns true if successful
+ * and nil otherwise. If wantsCoordinates is true, an array of arrays of the coordinates in the
+ * removed frames is returned if operation is successful.
*/
static VALUE
-s_Molecule_MeasureBond(VALUE self, VALUE nval1, VALUE nval2)
+s_Molecule_RemoveFrames(int argc, VALUE *argv, VALUE self)
{
+ VALUE val, flag;
+ VALUE retval;
Molecule *mol;
- Vector v1, v2;
+ IntGroup *ig;
+ int count;
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));
+ rb_scan_args(argc, argv, "11", &val, &flag);
+ ig = IntGroupFromValue(val);
+ count = IntGroupGetCount(ig);
+ if (RTEST(flag)) {
+ /* Create return value before removing frames */
+ VALUE coords;
+ int i, j, n;
+ Atom *ap;
+ Vector v;
+ retval = rb_ary_new2(count);
+ for (i = 0; i < count; i++) {
+ n = IntGroupGetNthPoint(ig, i);
+ coords = rb_ary_new2(mol->natoms);
+ for (j = 0, ap = mol->atoms; j < mol->natoms; j++, ap = ATOM_NEXT(ap)) {
+ if (n < ap->nframes && n != mol->cframe)
+ v = ap->frames[n];
+ else v = ap->r;
+ rb_ary_push(coords, ValueFromVector(&v));
+ }
+ rb_ary_push(retval, coords);
+ }
+ } else retval = Qtrue;
+ if (MolActionCreateAndPerform(mol, gMolActionRemoveFrames, ig) >= 0)
+ return retval;
+ else return Qnil;
}
/*
* call-seq:
- * measure_angle(n1, n2, n3) -> Float
+ * each_frame {|n| ...}
*
- * 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.
+ * Set the frame number from 0 to nframes-1 and execute the block. The block argument is
+ * the frame number. After completion, the original frame number is restored.
*/
static VALUE
-s_Molecule_MeasureAngle(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3)
+s_Molecule_EachFrame(VALUE self)
{
+ int i, cframe, nframes;
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);
+ cframe = mol->cframe;
+ nframes = MoleculeGetNumberOfFrames(mol);
+ if (nframes > 0) {
+ for (i = 0; i < nframes; i++) {
+ MoleculeSelectFrame(mol, i, 1);
+ rb_yield(INT2NUM(i));
+ }
+ MoleculeSelectFrame(mol, cframe, 1);
+ }
+ return self;
}
/*
* call-seq:
- * measure_dihedral(n1, n2, n3, n4) -> Float
+ * get_coord_from_frame(index, group = nil)
*
- * 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.
+ * Copy the coordinates from the indicated frame. If group is specified, only the specified atoms
+ * are modified. Third argument (cflag) is now obsolete (it used to specify whether the cell parameters are to be
+ * copied; now they are always copied)
*/
static VALUE
-s_Molecule_MeasureDihedral(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3, VALUE nval4)
+s_Molecule_GetCoordFromFrame(int argc, VALUE *argv, VALUE self)
{
- Molecule *mol;
- Vector v1, v2, v3, v4;
- Double d;
+ Molecule *mol;
+ VALUE ival, gval, cval;
+ Int index, i, j, n, nn;
+ IntGroup *ig;
+ IntGroupIterator iter;
+ Atom *ap;
+ Vector *vp;
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);
+ rb_scan_args(argc, argv, "12", &ival, &gval, &cval);
+ if (argc == 3)
+ rb_warn("The 3rd argument to get_coord_from_frame() is now obsolete");
+ index = NUM2INT(rb_Integer(ival));
+ if (index < 0 || index >= (n = MoleculeGetNumberOfFrames(mol))) {
+ if (n == 0)
+ rb_raise(rb_eMolbyError, "No frame is present");
+ else
+ rb_raise(rb_eMolbyError, "Frame index (%d) out of range (should be 0..%d)", index, n - 1);
+ }
+ if (gval == Qnil) {
+ ig = IntGroupNewWithPoints(0, mol->natoms, -1);
+ } else {
+ ig = s_Molecule_AtomGroupFromValue(self, gval);
+ }
+ n = IntGroupGetCount(ig);
+ if (n > 0) {
+ vp = (Vector *)calloc(sizeof(Vector), n);
+ IntGroupIteratorInit(ig, &iter);
+ j = 0;
+ nn = 0;
+ while ((i = IntGroupIteratorNext(&iter)) >= 0) {
+ ap = ATOM_AT_INDEX(mol->atoms, i);
+ if (index < ap->nframes) {
+ vp[j] = ap->frames[index];
+ nn++;
+ } else {
+ vp[j] = ap->r;
+ }
+ j++;
+ }
+ if (nn > 0)
+ MolActionCreateAndPerform(mol, gMolActionSetAtomPositions, ig, n, vp);
+ free(vp);
+ if (mol->cell != NULL && mol->frame_cells != NULL && index < mol->nframe_cells) {
+ vp = mol->frame_cells + index * 4;
+ MolActionCreateAndPerform(mol, gMolActionSetBox, vp, vp + 1, vp + 2, vp + 3, -1, 0);
+ }
+ IntGroupIteratorRelease(&iter);
+ }
+ /* Copy the extra properties */
+ IntGroupRelease(ig);
+ for (i = 0; i < mol->nmolprops; i++) {
+ Double *dp = (Double *)malloc(sizeof(Double));
+ ig = IntGroupNew();
+ IntGroupAdd(ig, mol->cframe, 1);
+ *dp = mol->molprops[i].propvals[index];
+ MolActionCreateAndPerform(mol, gMolActionSetProperty, i, ig, 1, dp);
+ free(dp);
+ IntGroupRelease(ig);
+ }
+
+ return self;
}
/*
* call-seq:
- * expand_by_symmetry(group, sym, dx=0, dy=0, dz=0, allow_overlap = false) -> Array
+ * reorder_frames(old_indices)
*
- * Expand the specified part of the molecule by the given symmetry operation.
- * Returns the array of atom indices corresponding to the expanded atoms.
- * If allow_overlap is true, then new atoms are created even when the
- * coordinates coincide with the some other atom (special position) of the
- * same element; otherwise, such atom will not be created and the index of the
- * existing atom is given in the returned array.
+ * Reorder the frames. The argument is an array of integers that specify the 'old'
+ * frame numbers. Thus, if the argument is [2,0,1], then the new frames 0/1/2 are the
+ * same as the old frames 2/0/1, respectively.
+ * The argument must have the same number of integers as the number of frames.
*/
static VALUE
-s_Molecule_ExpandBySymmetry(int argc, VALUE *argv, VALUE self)
+s_Molecule_ReorderFrames(VALUE self, VALUE aval)
{
- Molecule *mol;
- VALUE gval, sval, xval, yval, zval, rval, oval;
- IntGroup *ig;
- Int n[4], allow_overlap;
- Int natoms;
- Int nidx, *idx;
-
+ Molecule *mol;
+ Int *ip, *ip2, i, n, nframes;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "24", &gval, &sval, &xval, &yval, &zval, &oval);
- 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)));
- allow_overlap = (RTEST(oval) ? 1 : 0);
- 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], allow_overlap, &nidx, &idx);
-
- rval = rb_ary_new2(nidx);
- while (--nidx >= 0) {
- rb_ary_store(rval, nidx, INT2NUM(idx[nidx]));
+ aval = rb_ary_to_ary(aval);
+ nframes = MoleculeGetNumberOfFrames(mol);
+ if (RARRAY_LEN(aval) != nframes)
+ rb_raise(rb_eMolbyError, "The argument must have the same number of integers as the number of frames");
+ ip2 = (Int *)calloc(sizeof(Int), nframes);
+ ip = (Int *)calloc(sizeof(Int), nframes);
+ for (i = 0; i < nframes; i++) {
+ n = NUM2INT(rb_Integer(RARRAY_PTR(aval)[i]));
+ if (n < 0 || n >= nframes || ip2[n] != 0) {
+ free(ip2);
+ free(ip);
+ if (n < 0 || n >= nframes)
+ rb_raise(rb_eMolbyError, "The argument (%d) is out of range", n);
+ else
+ rb_raise(rb_eMolbyError, "The argument has duplicated entry (%d)", n);
+ }
+ ip2[n] = 1;
+ ip[i] = n;
}
-/* 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;
+ free(ip2);
+ MolActionCreateAndPerform(mol, gMolActionReorderFrames, nframes, ip);
+ free(ip);
+ return self;
}
+#pragma mark ------ Fragments ------
+
/*
* call-seq:
- * amend_by_symmetry(group = nil) -> IntGroup
+ * fragment(n1, *exatoms) -> IntGroup
+ * fragment(group, *exatoms) -> IntGroup
*
- * Expand the specified part of the molecule by the given symmetry operation.
- * Returns an IntGroup containing the added atoms.
+ * Get the fragment including the atom n1 or the atom group. If additional arguments are given,
+ * those atoms will not be counted during the search.
*/
static VALUE
-s_Molecule_AmendBySymmetry(int argc, VALUE *argv, VALUE self)
+s_Molecule_Fragment(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- IntGroup *ig, *ig2;
- VALUE rval, gval;
+ IntGroup *baseg, *ig, *exatoms;
+ int n;
+ volatile VALUE nval, exval;
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;
+ rb_scan_args(argc, argv, "1*", &nval, &exval);
+ if (rb_obj_is_kind_of(nval, rb_cNumeric) || rb_obj_is_kind_of(nval, rb_cString)) {
+ baseg = NULL;
+ n = NUM2INT(s_Molecule_AtomIndex(self, nval));
+ } else {
+ baseg = s_Molecule_AtomGroupFromValue(self, nval);
+ }
+ if (RARRAY_LEN(exval) == 0) {
+ exatoms = NULL;
+ } else {
+ exval = s_Molecule_AtomGroup(RARRAY_LEN(exval), RARRAY_PTR(exval), self);
+ Data_Get_Struct(exval, IntGroup, exatoms);
+ }
+ if (baseg == NULL) {
+ ig = MoleculeFragmentExcludingAtomGroup(mol, n, exatoms);
+ } else {
+ IntGroupIterator iter;
+ IntGroupIteratorInit(baseg, &iter);
+ if ((n = IntGroupIteratorNext(&iter)) < 0) {
+ ig = IntGroupNew();
+ } else {
+ ig = MoleculeFragmentExcludingAtomGroup(mol, n, exatoms);
+ if (ig != NULL) {
+ while ((n = IntGroupIteratorNext(&iter)) >= 0) {
+ IntGroup *subg;
+ subg = MoleculeFragmentExcludingAtomGroup(mol, n, exatoms);
+ if (subg != NULL) {
+ IntGroupAddIntGroup(ig, subg);
+ IntGroupRelease(subg);
+ }
+ }
+ }
+ }
+ IntGroupIteratorRelease(&iter);
+ IntGroupRelease(baseg);
+ }
+ if (ig == NULL)
+ rb_raise(rb_eMolbyError, "invalid specification of molecular fragment");
+ nval = ValueFromIntGroup(ig);
+ IntGroupRelease(ig);
+ return nval;
}
/*
* call-seq:
- * transform_for_symop(symop, is_cartesian = nil) -> Transform
+ * fragments(exclude = nil)
*
- * 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.
+ * Returns the fragments as an array of IntGroups.
+ * If exclude is given (as an array or an IntGroup), then those atoms are excluded
+ * in defining the fragment.
*/
static VALUE
-s_Molecule_TransformForSymop(int argc, VALUE *argv, VALUE self)
+s_Molecule_Fragments(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE sval, fval;
- Symop symop;
- Transform tr;
+ IntGroup *ag, *fg, *eg;
+ VALUE gval, exval, retval;
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 (mol == NULL)
+ return Qnil;
+ if (mol->natoms == 0)
+ return rb_ary_new();
+ rb_scan_args(argc, argv, "01", &exval);
+ if (exval == Qnil)
+ eg = NULL;
+ else
+ eg = IntGroupFromValue(exval);
+ ag = IntGroupNewWithPoints(0, mol->natoms, -1);
+ if (eg != NULL)
+ IntGroupRemoveIntGroup(ag, eg);
+ retval = rb_ary_new();
+ while (IntGroupGetCount(ag) > 0) {
+ int n = IntGroupGetNthPoint(ag, 0);
+ fg = MoleculeFragmentExcludingAtomGroup(mol, n, eg);
+ if (fg == NULL)
+ rb_raise(rb_eMolbyError, "internal error during each_fragment");
+ gval = ValueFromIntGroup(fg);
+ rb_ary_push(retval, gval);
+ IntGroupRemoveIntGroup(ag, fg);
+ IntGroupRelease(fg);
}
- 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);
+ IntGroupRelease(ag);
+ if (eg != NULL)
+ IntGroupRelease(eg);
+ return retval;
}
-
+
/*
* call-seq:
- * symop_for_transform(transform, is_cartesian = nil) -> [sym, dx, dy, dz]
+ * each_fragment(exclude = nil) {|group| ...}
*
- * 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.
+ * Execute the block, with the IntGroup object for each fragment as the argument.
+ * Atoms or bonds should not be added or removed during the execution of the block.
+ * If exclude is given (as an array or an IntGroup), then those atoms are excluded
+ * in defining the fragment.
*/
static VALUE
-s_Molecule_SymopForTransform(int argc, VALUE *argv, VALUE self)
+s_Molecule_EachFragment(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE tval, fval;
- Symop symop;
- Transform tr;
- int n;
+ IntGroup *ag, *fg, *eg;
+ VALUE gval, exval;
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 */
+ if (mol == NULL || mol->natoms == 0)
+ return self;
+ rb_scan_args(argc, argv, "01", &exval);
+ if (exval == Qnil)
+ eg = NULL;
+ else
+ eg = IntGroupFromValue(exval);
+ ag = IntGroupNewWithPoints(0, mol->natoms, -1);
+ if (eg != NULL)
+ IntGroupRemoveIntGroup(ag, eg);
+ while (IntGroupGetCount(ag) > 0) {
+ int n = IntGroupGetNthPoint(ag, 0);
+ fg = MoleculeFragmentExcludingAtomGroup(mol, n, eg);
+ if (fg == NULL)
+ rb_raise(rb_eMolbyError, "internal error during each_fragment");
+ gval = ValueFromIntGroup(fg);
+ rb_yield(gval);
+ IntGroupRemoveIntGroup(ag, fg);
+ IntGroupRelease(fg);
}
+ IntGroupRelease(ag);
+ if (eg != NULL)
+ IntGroupRelease(eg);
+ return self;
}
/*
* call-seq:
- * wrap_unit_cell(group) -> Vector3D
+ * detachable?(group) -> [n1, n2]
*
- * 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.
+ * Check whether the group is 'detachable', i.e. the group is bound to the rest
+ * of the molecule via only one bond. If it is, then the indices of the atoms
+ * belonging to the bond is returned, the first element being the atom included
+ * in the fragment. Otherwise, Qnil is returned.
*/
static VALUE
-s_Molecule_WrapUnitCell(VALUE self, VALUE gval)
+s_Molecule_Detachable_P(VALUE self, VALUE gval)
{
- Molecule *mol;
+ Molecule *mol;
IntGroup *ig;
- Vector v, cv, dv;
+ int n1, n2;
+ VALUE retval;
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);
+ if (MoleculeIsFragmentDetachable(mol, ig, &n1, &n2)) {
+ retval = rb_ary_new3(2, INT2NUM(n1), INT2NUM(n2));
+ } else retval = Qnil;
IntGroupRelease(ig);
- return ValueFromVector(&dv);
+ return retval;
}
/*
* call-seq:
- * find_conflicts(limit[, group1[, group2 [, ignore_exclusion]]]) -> [[n1, n2], [n3, n4], ...]
+ * bonds_on_border(group = selection) -> Array of Array of two Integers
*
- * 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.
+ * 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_FindConflicts(int argc, VALUE *argv, VALUE self)
+s_Molecule_BondsOnBorder(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;
-
+ Molecule *mol;
+ IntGroup *ig, *bg;
+ VALUE gval, retval;
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;
+ rb_scan_args(argc, argv, "01", &gval);
+ if (gval == Qnil) {
+ ig = MoleculeGetSelection(mol);
+ if (ig != NULL)
+ IntGroupRetain(ig);
} else {
- exinfo = NULL;
- exlist = NULL;
+ ig = s_Molecule_AtomGroupFromValue(self, gval);
}
- 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);
- }
+ 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);
}
- IntGroupIteratorRelease(&iter2);
- IntGroupIteratorRelease(&iter1);
- IntGroupRelease(ig2);
- IntGroupRelease(ig1);
- rval = rb_ary_new2(npairs);
- if (pairs != NULL) {
- for (i = 0; i < npairs; i++) {
- rb_ary_push(rval, rb_ary_new3(2, INT2NUM(pairs[i * 2]), INT2NUM(pairs[i * 2 + 1])));
- }
- free(pairs);
- }
- return rval;
+ IntGroupRelease(bg);
+ IntGroupRelease(ig);
+ return retval;
}
/* Calculate the transform that moves the current coordinates to the reference
return Qnil; /* Not reached */
}
+#pragma mark ------ Screen Display ------
+
/*
* call-seq:
* display
/*
* call-seq:
+ * hidden_atoms -> IntGroup
+ *
+ * Returns the currently hidden atoms.
+ */
+static VALUE
+s_Molecule_HiddenAtoms(VALUE self)
+{
+ rb_raise(rb_eMolbyError, "set_hidden_atoms is now obsolete. Try using Molecule#is_atom_visible or AtomRef#hidden.");
+ return Qnil; /* Not reached */
+}
+
+/*
+ * call-seq:
+ * set_hidden_atoms(IntGroup)
+ * self.hidden_atoms = IntGroup
+ *
+ * Hide the specified atoms. This operation is _not_ undoable.
+ */
+static VALUE
+s_Molecule_SetHiddenAtoms(VALUE self, VALUE val)
+{
+ rb_raise(rb_eMolbyError, "set_hidden_atoms is now obsolete. Try using Molecule#is_atom_visible or AtomRef#hidden.");
+ return Qnil; /* Not reached */
+}
+
+/*
+ * call-seq:
* show_graphite -> Integer
* show_graphite = Integer
* show_graphite = boolean
/*
* call-seq:
+ * export_graphic(fname, scale = 1.0, bg_color = -1)
+ *
+ * Export the current graphic to a PNG or TIF file (determined by the extension).
+ * bg_color: -1, same as screen; 0, transparent; 1, black; 2, white.
+ *
+ */
+static VALUE
+s_Molecule_ExportGraphic(int argc, VALUE *argv, VALUE self)
+{
+ Molecule *mol;
+ VALUE fval, sval, bval;
+ char *fname;
+ float scale;
+ int bg_color;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->mview == NULL)
+ rb_raise(rb_eMolbyError, "The molecule has no associated graphic view");
+ rb_scan_args(argc, argv, "12", &fval, &sval, &bval);
+ fname = FileStringValuePtr(fval);
+ if (sval == Qnil)
+ scale = 1.0;
+ else scale = NUM2DBL(rb_Float(sval));
+ if (bval == Qnil)
+ bg_color = -1;
+ else bg_color = NUM2INT(rb_Integer(bval));
+ if (MainViewCallback_exportGraphic(mol->mview, fname, scale, bg_color) == 0)
+ return fval;
+ else return Qnil;
+}
+
+#pragma mark ------ Graphics ------
+
+/*
+ * call-seq:
* insert_graphic(index, kind, color, points, fill = nil) -> integer
*
* Create a new graphic object and insert at the given graphic index (if -1, then append at the last).
return Qnil;
}
+#pragma mark ------ MD Support ------
+
/*
* call-seq:
* md_arena -> MDArena
/*
* call-seq:
+ * start_step -> Integer
+ *
+ * Returns the start step (defined by dcd format).
+ */
+static VALUE
+s_Molecule_StartStep(VALUE self)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ return INT2NUM(mol->startStep);
+}
+
+/*
+ * call-seq:
+ * start_step = Integer
+ *
+ * Set the start step (defined by dcd format).
+ */
+static VALUE
+s_Molecule_SetStartStep(VALUE self, VALUE val)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ mol->startStep = NUM2INT(rb_Integer(val));
+ return val;
+}
+
+/*
+ * call-seq:
+ * steps_per_frame -> Integer
+ *
+ * Returns the number of steps between frames (defined by dcd format).
+ */
+static VALUE
+s_Molecule_StepsPerFrame(VALUE self)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ return INT2NUM(mol->stepsPerFrame);
+}
+
+/*
+ * call-seq:
+ * steps_per_frame = Integer
+ *
+ * Set the number of steps between frames (defined by dcd format).
+ */
+static VALUE
+s_Molecule_SetStepsPerFrame(VALUE self, VALUE val)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ mol->stepsPerFrame = NUM2INT(rb_Integer(val));
+ return val;
+}
+
+/*
+ * call-seq:
+ * ps_per_step -> Float
+ *
+ * Returns the time increment (in picoseconds) for one step (defined by dcd format).
+ */
+static VALUE
+s_Molecule_PsPerStep(VALUE self)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ return rb_float_new(mol->psPerStep);
+}
+
+/*
+ * call-seq:
+ * ps_per_step = Float
+ *
+ * Set the time increment (in picoseconds) for one step (defined by dcd format).
+ */
+static VALUE
+s_Molecule_SetPsPerStep(VALUE self, VALUE val)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ mol->psPerStep = NUM2DBL(rb_Float(val));
+ return val;
+}
+
+static VALUE
+s_Molecule_BondParIsObsolete(VALUE self, VALUE val)
+{
+ rb_raise(rb_eMolbyError, "Molecule#bond_par, angle_par, dihedral_par, improper_par, vdw_par are now obsolete. You can use MDArena#bond_par, angle_par, dihedral_par, improper_par, vdw_par instead, and probably these are what you really want.");
+}
+
+#pragma mark ------ MO Handling ------
+
+/*
+ * call-seq:
* selectedMO -> IntGroup
*
* Returns a group of selected mo in the "MO Info" table. If the MO info table
*
* Returns either "RHF", "UHF", or "ROHF". If no MO info is present, returns nil.
*/
-static VALUE
-s_Molecule_MOType(VALUE self)
-{
- return s_Molecule_GetMOInfo(self, sTypeSym);
-}
-
-
-#if 0
-/*
- * call-seq:
- * allocate_basis_set_record(rflag, ne_alpha, ne_beta)
- *
- * To be used internally. Allocate a basis set record. rflag: 0, unrestricted; 1, restricted.
- * ne_alpha, ne_beta: number of alpha/beta electrons.
- */
-static VALUE
-s_Molecule_AllocateBasisSetRecord(VALUE self, VALUE rval, VALUE naval, VALUE nbval)
-{
- Molecule *mol;
- Int rflag, na, nb, n;
- Data_Get_Struct(self, Molecule, mol);
- rflag = NUM2INT(rb_Integer(rval));
- na = NUM2INT(rb_Integer(naval));
- nb = NUM2INT(rb_Integer(nbval));
- n = MoleculeSetMOInfo(mol, rflag, na, nb);
- if (n == -1)
- rb_raise(rb_eMolbyError, "Molecule is emptry");
- else if (n == -2)
- rb_raise(rb_eMolbyError, "Low memory");
- else if (n != 0)
- rb_raise(rb_eMolbyError, "Unknown error");
- return self;
+static VALUE
+s_Molecule_MOType(VALUE self)
+{
+ return s_Molecule_GetMOInfo(self, sTypeSym);
}
-#endif
+
+#pragma mark ------ Molecular Topology ------
/*
* call-seq:
return Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
}
+#pragma mark ------ Molecular Properties ------
+
/*
* call-seq:
* set_property(name, value[, index]) -> value
return rval;
}
-/*
- * call-seq:
- * export_graphic(fname, scale = 1.0, bg_color = -1)
- *
- * Export the current graphic to a PNG or TIF file (determined by the extension).
- * bg_color: -1, same as screen; 0, transparent; 1, black; 2, white.
- *
- */
-static VALUE
-s_Molecule_ExportGraphic(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE fval, sval, bval;
- char *fname;
- float scale;
- int bg_color;
- Data_Get_Struct(self, Molecule, mol);
- if (mol->mview == NULL)
- rb_raise(rb_eMolbyError, "The molecule has no associated graphic view");
- rb_scan_args(argc, argv, "12", &fval, &sval, &bval);
- fname = FileStringValuePtr(fval);
- if (sval == Qnil)
- scale = 1.0;
- else scale = NUM2DBL(rb_Float(sval));
- if (bval == Qnil)
- bg_color = -1;
- else bg_color = NUM2INT(rb_Integer(bval));
- if (MainViewCallback_exportGraphic(mol->mview, fname, scale, bg_color) == 0)
- return fval;
- else return Qnil;
-}
+#pragma mark ------ Class methods ------
/*
* call-seq:
return ary;
}
-/*
- * call-seq:
- * error_message -> String
- *
- * Get the error_message from the last load/save method. If no error, returns nil.
- */
-static VALUE
-s_Molecule_ErrorMessage(VALUE klass)
-{
- if (gLoadSaveErrorMessage == NULL)
- return Qnil;
- else return rb_str_new2(gLoadSaveErrorMessage);
-}
-
-/*
- * call-seq:
- * set_error_message(String)
- * Molecule.error_message = String
- *
- * Get the error_message from the last load/save method. If no error, returns nil.
- */
-static VALUE
-s_Molecule_SetErrorMessage(VALUE klass, VALUE sval)
-{
- if (gLoadSaveErrorMessage != NULL) {
- free(gLoadSaveErrorMessage);
- gLoadSaveErrorMessage = NULL;
- }
- if (sval != Qnil) {
- sval = rb_str_to_str(sval);
- gLoadSaveErrorMessage = strdup(StringValuePtr(sval));
- }
- return sval;
-}
-
-/*
- * call-seq:
- * self == Molecule -> boolean
- *
- * True if the two arguments point to the same molecule.
- */
-static VALUE
-s_Molecule_Equal(VALUE self, VALUE val)
-{
- if (rb_obj_is_kind_of(val, rb_cMolecule)) {
- Molecule *mol1, *mol2;
- Data_Get_Struct(self, Molecule, mol1);
- Data_Get_Struct(val, Molecule, mol2);
- return (mol1 == mol2 ? Qtrue : Qfalse);
- } else return Qfalse;
-}
+#pragma mark ------ Call Subprocess ------
/* The callback functions for call_subprocess_async */
static int
return INT2NUM(n);
}
+#pragma mark ====== Define Molby Classes ======
+
void
Init_Molby(void)
{
/* class Molecule */
rb_cMolecule = rb_define_class_under(rb_mMolby, "Molecule", rb_cObject);
+
rb_define_alloc_func(rb_cMolecule, s_Molecule_Alloc);
rb_define_private_method(rb_cMolecule, "initialize", s_Molecule_Initialize, -1);
rb_define_private_method(rb_cMolecule, "initialize_copy", s_Molecule_InitCopy, 1);
+ rb_define_method(rb_cMolecule, "atom_index", s_Molecule_AtomIndex, 1);
+ rb_define_method(rb_cMolecule, "==", s_Molecule_Equal, 1);
+
rb_define_method(rb_cMolecule, "loadmbsf", s_Molecule_Loadmbsf, -1);
rb_define_method(rb_cMolecule, "loadpsf", s_Molecule_Loadpsf, -1);
rb_define_alias(rb_cMolecule, "loadpsfx", "loadpsf");
rb_define_method(rb_cMolecule, "loadfchk", s_Molecule_Loadfchk, -1);
rb_define_alias(rb_cMolecule, "loadfch", "loadfchk");
rb_define_method(rb_cMolecule, "loaddat", s_Molecule_Loaddat, -1);
- rb_define_method(rb_cMolecule, "molload", s_Molecule_Load, -1);
- rb_define_method(rb_cMolecule, "molsave", s_Molecule_Save, -1);
rb_define_method(rb_cMolecule, "savembsf", s_Molecule_Savembsf, 1);
rb_define_method(rb_cMolecule, "savepsf", s_Molecule_Savepsf, 1);
rb_define_alias(rb_cMolecule, "savepsfx", "savepsf");
rb_define_method(rb_cMolecule, "savepdb", s_Molecule_Savepdb, 1);
rb_define_method(rb_cMolecule, "savedcd", s_Molecule_Savedcd, 1);
-/* rb_define_method(rb_cMolecule, "savetep", s_Molecule_Savetep, 1); */
+ rb_define_method(rb_cMolecule, "molload", s_Molecule_Load, -1);
+ rb_define_method(rb_cMolecule, "molsave", s_Molecule_Save, -1);
+ rb_define_singleton_method(rb_cMolecule, "open", s_Molecule_Open, -1);
+ rb_define_singleton_method(rb_cMolecule, "error_message", s_Molecule_ErrorMessage, 0);
+ rb_define_singleton_method(rb_cMolecule, "set_error_message", s_Molecule_SetErrorMessage, 1);
+ rb_define_singleton_method(rb_cMolecule, "error_message=", s_Molecule_SetErrorMessage, 1);
+
rb_define_method(rb_cMolecule, "name", s_Molecule_Name, 0);
rb_define_method(rb_cMolecule, "set_name", s_Molecule_SetName, 1);
rb_define_method(rb_cMolecule, "path", s_Molecule_Path, 0);
rb_define_method(rb_cMolecule, "dir", s_Molecule_Dir, 0);
rb_define_method(rb_cMolecule, "inspect", s_Molecule_Inspect, 0);
+
rb_define_method(rb_cMolecule, "atoms", s_Molecule_Atoms, 0);
rb_define_method(rb_cMolecule, "bonds", s_Molecule_Bonds, 0);
rb_define_method(rb_cMolecule, "angles", s_Molecule_Angles, 0);
rb_define_method(rb_cMolecule, "nangles", s_Molecule_Nangles, 0);
rb_define_method(rb_cMolecule, "ndihedrals", s_Molecule_Ndihedrals, 0);
rb_define_method(rb_cMolecule, "nimpropers", s_Molecule_Nimpropers, 0);
-
- rb_define_method(rb_cMolecule, "bond_par", s_Molecule_BondParIsObsolete, 1);
- rb_define_method(rb_cMolecule, "angle_par", s_Molecule_BondParIsObsolete, 1);
- rb_define_method(rb_cMolecule, "dihedral_par", s_Molecule_BondParIsObsolete, 1);
- rb_define_method(rb_cMolecule, "improper_par", s_Molecule_BondParIsObsolete, 1);
- rb_define_method(rb_cMolecule, "vdw_par", s_Molecule_BondParIsObsolete, 1);
-
- rb_define_method(rb_cMolecule, "start_step", s_Molecule_StartStep, 0);
- rb_define_method(rb_cMolecule, "start_step=", s_Molecule_SetStartStep, 1);
- rb_define_method(rb_cMolecule, "steps_per_frame", s_Molecule_StepsPerFrame, 0);
- rb_define_method(rb_cMolecule, "steps_per_frame=", s_Molecule_SetStepsPerFrame, 1);
- 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, "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, "min_residue_number", s_Molecule_MinResSeq, -1);
+
rb_define_method(rb_cMolecule, "each_atom", s_Molecule_EachAtom, -1);
- rb_define_method(rb_cMolecule, "cell", s_Molecule_Cell, 0);
- rb_define_method(rb_cMolecule, "cell=", s_Molecule_SetCell, -1);
- rb_define_alias(rb_cMolecule, "set_cell", "cell=");
- rb_define_method(rb_cMolecule, "box", s_Molecule_Box, 0);
- rb_define_method(rb_cMolecule, "box=", s_Molecule_SetBox, 1);
- rb_define_method(rb_cMolecule, "set_box", s_Molecule_SetBox, -2);
- rb_define_method(rb_cMolecule, "cell_transform", s_Molecule_CellTransform, 0);
- rb_define_method(rb_cMolecule, "cell_periodicity", s_Molecule_CellPeriodicity, 0);
- rb_define_method(rb_cMolecule, "cell_periodicity=", s_Molecule_SetCellPeriodicity, 1);
- rb_define_alias(rb_cMolecule, "set_cell_periodicity", "cell_periodicity=");
- rb_define_method(rb_cMolecule, "cell_flexibility", s_Molecule_CellFlexibility, 0);
- rb_define_method(rb_cMolecule, "cell_flexibility=", s_Molecule_SetCellFlexibility, 1);
- rb_define_alias(rb_cMolecule, "set_cell_flexibility", "cell_flexibility=");
- rb_define_method(rb_cMolecule, "symmetry", s_Molecule_Symmetry, 0);
- rb_define_alias(rb_cMolecule, "symmetries", "symmetry");
- rb_define_method(rb_cMolecule, "nsymmetries", s_Molecule_Nsymmetries, 0);
- rb_define_method(rb_cMolecule, "add_symmetry", s_Molecule_AddSymmetry, 1);
- rb_define_method(rb_cMolecule, "remove_symmetry", s_Molecule_RemoveSymmetry, -1);
- rb_define_alias(rb_cMolecule, "remove_symmetries", "remove_symmetry");
+ rb_define_method(rb_cMolecule, "atom_group", s_Molecule_AtomGroup, -1);
+ rb_define_method(rb_cMolecule, "selection", s_Molecule_Selection, 0);
+ rb_define_method(rb_cMolecule, "selection=", s_Molecule_SetSelection, 1);
+ rb_define_method(rb_cMolecule, "set_undoable_selection", s_Molecule_SetUndoableSelection, 1);
+
rb_define_method(rb_cMolecule, "extract", s_Molecule_Extract, -1);
rb_define_method(rb_cMolecule, "add", s_Molecule_Add, 1);
rb_define_alias(rb_cMolecule, "+", "add");
rb_define_method(rb_cMolecule, "remove", s_Molecule_Remove, 1);
- rb_define_method(rb_cMolecule, "atom_group", s_Molecule_AtomGroup, -1);
- rb_define_method(rb_cMolecule, "atom_index", s_Molecule_AtomIndex, 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_method(rb_cMolecule, "assign_residue", s_Molecule_AssignResidue, 2);
rb_define_method(rb_cMolecule, "offset_residue", s_Molecule_OffsetResidue, 2);
rb_define_method(rb_cMolecule, "renumber_atoms", s_Molecule_RenumberAtoms, 1);
+
+ rb_define_method(rb_cMolecule, "set_atom_attr", s_Molecule_SetAtomAttr, 3);
+ rb_define_method(rb_cMolecule, "get_atom_attr", s_Molecule_GetAtomAttr, 2);
+ rb_define_method(rb_cMolecule, "register_undo", s_Molecule_RegisterUndo, -1);
+ rb_define_method(rb_cMolecule, "undo_enabled?", s_Molecule_UndoEnabled, 0);
+ rb_define_method(rb_cMolecule, "undo_enabled=", s_Molecule_SetUndoEnabled, 1);
+
+ rb_define_method(rb_cMolecule, "center_of_mass", s_Molecule_CenterOfMass, -1);
+ rb_define_method(rb_cMolecule, "centralize", s_Molecule_Centralize, -1);
+ rb_define_method(rb_cMolecule, "bounds", s_Molecule_Bounds, -1);
+ rb_define_method(rb_cMolecule, "measure_bond", s_Molecule_MeasureBond, 2);
+ rb_define_method(rb_cMolecule, "measure_angle", s_Molecule_MeasureAngle, 3);
+ rb_define_method(rb_cMolecule, "measure_dihedral", s_Molecule_MeasureDihedral, 4);
+ rb_define_method(rb_cMolecule, "find_conflicts", s_Molecule_FindConflicts, -1);
rb_define_method(rb_cMolecule, "find_close_atoms", s_Molecule_FindCloseAtoms, -1);
rb_define_method(rb_cMolecule, "guess_bonds", s_Molecule_GuessBonds, -1);
- rb_define_method(rb_cMolecule, "selection", s_Molecule_Selection, 0);
- rb_define_method(rb_cMolecule, "selection=", s_Molecule_SetSelection, 1);
- rb_define_method(rb_cMolecule, "set_undoable_selection", s_Molecule_SetUndoableSelection, 1);
- rb_define_method(rb_cMolecule, "hidden_atoms", s_Molecule_HiddenAtoms, 0); /* obsolete */
- rb_define_method(rb_cMolecule, "hidden_atoms=", s_Molecule_SetHiddenAtoms, 1); /* obsolete */
- rb_define_alias(rb_cMolecule, "set_hidden_atoms", "hidden_atoms="); /* obsolete */
- rb_define_method(rb_cMolecule, "frame", s_Molecule_Frame, 0);
+
+ rb_define_method(rb_cMolecule, "cell", s_Molecule_Cell, 0);
+ rb_define_method(rb_cMolecule, "cell=", s_Molecule_SetCell, -1);
+ rb_define_alias(rb_cMolecule, "set_cell", "cell=");
+ rb_define_method(rb_cMolecule, "box", s_Molecule_Box, 0);
+ rb_define_method(rb_cMolecule, "box=", s_Molecule_SetBox, 1);
+ rb_define_method(rb_cMolecule, "set_box", s_Molecule_SetBox, -2);
+ rb_define_method(rb_cMolecule, "cell_periodicity", s_Molecule_CellPeriodicity, 0);
+ rb_define_method(rb_cMolecule, "cell_periodicity=", s_Molecule_SetCellPeriodicity, 1);
+ rb_define_alias(rb_cMolecule, "set_cell_periodicity", "cell_periodicity=");
+ rb_define_method(rb_cMolecule, "cell_flexibility", s_Molecule_CellFlexibility, 0);
+ rb_define_method(rb_cMolecule, "cell_flexibility=", s_Molecule_SetCellFlexibility, 1);
+ rb_define_alias(rb_cMolecule, "set_cell_flexibility", "cell_flexibility=");
+ rb_define_method(rb_cMolecule, "cell_transform", s_Molecule_CellTransform, 0);
+ rb_define_method(rb_cMolecule, "symmetry", s_Molecule_Symmetry, 0);
+ rb_define_alias(rb_cMolecule, "symmetries", "symmetry");
+ rb_define_method(rb_cMolecule, "nsymmetries", s_Molecule_Nsymmetries, 0);
+ rb_define_method(rb_cMolecule, "add_symmetry", s_Molecule_AddSymmetry, 1);
+ rb_define_method(rb_cMolecule, "remove_symmetry", s_Molecule_RemoveSymmetry, -1);
+ rb_define_alias(rb_cMolecule, "remove_symmetries", "remove_symmetry");
+ rb_define_method(rb_cMolecule, "wrap_unit_cell", s_Molecule_WrapUnitCell, 1);
+ rb_define_method(rb_cMolecule, "expand_by_symmetry", s_Molecule_ExpandBySymmetry, -1);
+ rb_define_method(rb_cMolecule, "amend_by_symmetry", s_Molecule_AmendBySymmetry, -1);
+
+ rb_define_method(rb_cMolecule, "translate", s_Molecule_Translate, -1);
+ rb_define_method(rb_cMolecule, "rotate", s_Molecule_Rotate, -1);
+ rb_define_method(rb_cMolecule, "reflect", s_Molecule_Reflect, -1);
+ rb_define_method(rb_cMolecule, "invert", s_Molecule_Invert, -1);
+ rb_define_method(rb_cMolecule, "transform", s_Molecule_Transform, -1);
+ rb_define_method(rb_cMolecule, "transform_for_symop", s_Molecule_TransformForSymop, -1);
+ rb_define_method(rb_cMolecule, "symop_for_transform", s_Molecule_SymopForTransform, -1);
+
rb_define_method(rb_cMolecule, "frame=", s_Molecule_SelectFrame, 1);
rb_define_alias(rb_cMolecule, "select_frame", "frame=");
+ rb_define_method(rb_cMolecule, "frame", s_Molecule_Frame, 0);
rb_define_method(rb_cMolecule, "nframes", s_Molecule_Nframes, 0);
- rb_define_method(rb_cMolecule, "create_frame", s_Molecule_CreateFrames, -1);
rb_define_method(rb_cMolecule, "insert_frame", s_Molecule_InsertFrames, -1);
+ rb_define_method(rb_cMolecule, "create_frame", s_Molecule_CreateFrames, -1);
rb_define_method(rb_cMolecule, "remove_frame", s_Molecule_RemoveFrames, -1);
rb_define_alias(rb_cMolecule, "create_frames", "create_frame");
rb_define_alias(rb_cMolecule, "insert_frames", "insert_frame");
rb_define_alias(rb_cMolecule, "remove_frames", "remove_frame");
- rb_define_method(rb_cMolecule, "reorder_frames", s_Molecule_ReorderFrames, 1);
rb_define_method(rb_cMolecule, "each_frame", s_Molecule_EachFrame, 0);
rb_define_method(rb_cMolecule, "get_coord_from_frame", s_Molecule_GetCoordFromFrame, -1);
- rb_define_method(rb_cMolecule, "register_undo", s_Molecule_RegisterUndo, -1);
- rb_define_method(rb_cMolecule, "undo_enabled?", s_Molecule_UndoEnabled, 0);
- rb_define_method(rb_cMolecule, "undo_enabled=", s_Molecule_SetUndoEnabled, 1);
- rb_define_method(rb_cMolecule, "set_atom_attr", s_Molecule_SetAtomAttr, 3);
- rb_define_method(rb_cMolecule, "get_atom_attr", s_Molecule_GetAtomAttr, 2);
+ rb_define_method(rb_cMolecule, "reorder_frames", s_Molecule_ReorderFrames, 1);
+
rb_define_method(rb_cMolecule, "fragment", s_Molecule_Fragment, -1);
rb_define_method(rb_cMolecule, "fragments", s_Molecule_Fragments, -1);
rb_define_method(rb_cMolecule, "each_fragment", s_Molecule_EachFragment, -1);
rb_define_method(rb_cMolecule, "detachable?", s_Molecule_Detachable_P, 1);
rb_define_method(rb_cMolecule, "bonds_on_border", s_Molecule_BondsOnBorder, -1);
- rb_define_method(rb_cMolecule, "translate", s_Molecule_Translate, -1);
- rb_define_method(rb_cMolecule, "rotate", s_Molecule_Rotate, -1);
- rb_define_method(rb_cMolecule, "reflect", s_Molecule_Reflect, -1);
- rb_define_method(rb_cMolecule, "invert", s_Molecule_Invert, -1);
- rb_define_method(rb_cMolecule, "transform", s_Molecule_Transform, -1);
- rb_define_method(rb_cMolecule, "center_of_mass", s_Molecule_CenterOfMass, -1);
- rb_define_method(rb_cMolecule, "centralize", s_Molecule_Centralize, -1);
- rb_define_method(rb_cMolecule, "bounds", s_Molecule_Bounds, -1);
- rb_define_method(rb_cMolecule, "measure_bond", s_Molecule_MeasureBond, 2);
- rb_define_method(rb_cMolecule, "measure_angle", s_Molecule_MeasureAngle, 3);
- rb_define_method(rb_cMolecule, "measure_dihedral", s_Molecule_MeasureDihedral, 4);
- rb_define_method(rb_cMolecule, "expand_by_symmetry", s_Molecule_ExpandBySymmetry, -1);
- rb_define_method(rb_cMolecule, "amend_by_symmetry", s_Molecule_AmendBySymmetry, -1);
- rb_define_method(rb_cMolecule, "transform_for_symop", s_Molecule_TransformForSymop, -1);
- rb_define_method(rb_cMolecule, "symop_for_transform", s_Molecule_SymopForTransform, -1);
- rb_define_method(rb_cMolecule, "wrap_unit_cell", s_Molecule_WrapUnitCell, 1);
- rb_define_method(rb_cMolecule, "find_conflicts", s_Molecule_FindConflicts, -1);
rb_define_method(rb_cMolecule, "fit_coordinates", s_Molecule_FitCoordinates, -1);
+
rb_define_method(rb_cMolecule, "display", s_Molecule_Display, 0);
rb_define_method(rb_cMolecule, "make_front", s_Molecule_MakeFront, 0);
rb_define_method(rb_cMolecule, "update_enabled?", s_Molecule_UpdateEnabled, 0);
rb_define_method(rb_cMolecule, "show_dummy_atoms=", s_Molecule_ShowDummyAtoms, -1);
rb_define_method(rb_cMolecule, "show_expanded", s_Molecule_ShowExpanded, -1);
rb_define_method(rb_cMolecule, "show_expanded=", s_Molecule_ShowExpanded, -1);
- rb_define_method(rb_cMolecule, "is_atom_visible", s_Molecule_IsAtomVisible, 1);
rb_define_method(rb_cMolecule, "show_ellipsoids", s_Molecule_ShowEllipsoids, -1);
rb_define_method(rb_cMolecule, "show_ellipsoids=", s_Molecule_ShowEllipsoids, -1);
+ rb_define_method(rb_cMolecule, "is_atom_visible", s_Molecule_IsAtomVisible, 1);
+ rb_define_method(rb_cMolecule, "hidden_atoms", s_Molecule_HiddenAtoms, 0); /* obsolete */
+ rb_define_method(rb_cMolecule, "hidden_atoms=", s_Molecule_SetHiddenAtoms, 1); /* obsolete */
+ rb_define_alias(rb_cMolecule, "set_hidden_atoms", "hidden_atoms="); /* obsolete */
rb_define_method(rb_cMolecule, "show_graphite", s_Molecule_ShowGraphite, -1);
rb_define_method(rb_cMolecule, "show_graphite=", s_Molecule_ShowGraphite, -1);
rb_define_method(rb_cMolecule, "show_graphite?", s_Molecule_ShowGraphiteFlag, 0);
rb_define_method(rb_cMolecule, "set_view_scale", s_Molecule_SetViewScale, 1);
rb_define_method(rb_cMolecule, "set_view_center", s_Molecule_SetViewCenter, 1);
rb_define_method(rb_cMolecule, "set_background_color", s_Molecule_SetBackgroundColor, -1);
- rb_define_method(rb_cMolecule, "create_graphic", s_Molecule_CreateGraphic, -1);
+ rb_define_method(rb_cMolecule, "export_graphic", s_Molecule_ExportGraphic, -1);
+
rb_define_method(rb_cMolecule, "insert_graphic", s_Molecule_InsertGraphic, -1);
+ rb_define_method(rb_cMolecule, "create_graphic", s_Molecule_CreateGraphic, -1);
rb_define_method(rb_cMolecule, "remove_graphic", s_Molecule_RemoveGraphic, 1);
rb_define_method(rb_cMolecule, "ngraphics", s_Molecule_NGraphics, 0);
rb_define_method(rb_cMolecule, "set_graphic_point", s_Molecule_SetGraphicPoint, 3);
rb_define_method(rb_cMolecule, "show_graphic", s_Molecule_ShowGraphic, 1);
rb_define_method(rb_cMolecule, "hide_graphic", s_Molecule_HideGraphic, 1);
rb_define_method(rb_cMolecule, "show_text", s_Molecule_ShowText, 1);
+
rb_define_method(rb_cMolecule, "md_arena", s_Molecule_MDArena, 0);
rb_define_method(rb_cMolecule, "set_parameter_attr", s_Molecule_SetParameterAttr, 5);
rb_define_method(rb_cMolecule, "parameter", s_Molecule_Parameter, 0);
+ rb_define_method(rb_cMolecule, "start_step", s_Molecule_StartStep, 0);
+ rb_define_method(rb_cMolecule, "start_step=", s_Molecule_SetStartStep, 1);
+ rb_define_method(rb_cMolecule, "steps_per_frame", s_Molecule_StepsPerFrame, 0);
+ rb_define_method(rb_cMolecule, "steps_per_frame=", s_Molecule_SetStepsPerFrame, 1);
+ 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, "bond_par", s_Molecule_BondParIsObsolete, 1); /* Obsolete */
+ rb_define_method(rb_cMolecule, "angle_par", s_Molecule_BondParIsObsolete, 1); /* Obsolete */
+ rb_define_method(rb_cMolecule, "dihedral_par", s_Molecule_BondParIsObsolete, 1); /* Obsolete */
+ rb_define_method(rb_cMolecule, "improper_par", s_Molecule_BondParIsObsolete, 1); /* Obsolete */
+ rb_define_method(rb_cMolecule, "vdw_par", s_Molecule_BondParIsObsolete, 1); /* Obsolete */
+
rb_define_method(rb_cMolecule, "selected_MO", s_Molecule_SelectedMO, 0);
rb_define_method(rb_cMolecule, "default_MO_grid", s_Molecule_GetDefaultMOGrid, -1);
rb_define_method(rb_cMolecule, "cubegen", s_Molecule_Cubegen, -1);
rb_define_method(rb_cMolecule, "get_gaussian_shell_info", s_Molecule_GetGaussianShellInfo, 1);
rb_define_method(rb_cMolecule, "get_gaussian_primitive_coefficients", s_Molecule_GetGaussianPrimitiveCoefficients, 1);
rb_define_method(rb_cMolecule, "get_gaussian_component_info", s_Molecule_GetGaussianComponentInfo, 1);
- rb_define_method(rb_cMolecule, "mo_type", s_Molecule_MOType, 0);
rb_define_method(rb_cMolecule, "clear_mo_coefficients", s_Molecule_ClearMOCoefficients, 0);
rb_define_method(rb_cMolecule, "set_mo_coefficients", s_Molecule_SetMOCoefficients, 3);
rb_define_method(rb_cMolecule, "get_mo_coefficients", s_Molecule_GetMOCoefficients, 1);
rb_define_method(rb_cMolecule, "get_mo_energy", s_Molecule_GetMOEnergy, 1);
- rb_define_method(rb_cMolecule, "get_mo_info", s_Molecule_GetMOInfo, 1);
rb_define_method(rb_cMolecule, "set_mo_info", s_Molecule_SetMOInfo, 1);
-/* rb_define_method(rb_cMolecule, "allocate_basis_set_record", s_Molecule_AllocateBasisSetRecord, 3); */
+ rb_define_method(rb_cMolecule, "get_mo_info", s_Molecule_GetMOInfo, 1);
+ rb_define_method(rb_cMolecule, "mo_type", s_Molecule_MOType, 0);
+
rb_define_method(rb_cMolecule, "search_equivalent_atoms", s_Molecule_SearchEquivalentAtoms, -1);
-
rb_define_method(rb_cMolecule, "create_pi_anchor", s_Molecule_CreatePiAnchor, -1);
+
rb_define_method(rb_cMolecule, "set_property", s_Molecule_SetProperty, -1);
rb_define_method(rb_cMolecule, "get_property", s_Molecule_GetProperty, -1);
rb_define_method(rb_cMolecule, "property_names", s_Molecule_PropertyNames, 0);
- rb_define_method(rb_cMolecule, "export_graphic", s_Molecule_ExportGraphic, -1);
- rb_define_method(rb_cMolecule, "==", s_Molecule_Equal, 1);
- rb_define_method(rb_cMolecule, "call_subprocess_async", s_Molecule_CallSubProcessAsync, -1);
-
rb_define_singleton_method(rb_cMolecule, "current", s_Molecule_Current, 0);
rb_define_singleton_method(rb_cMolecule, "[]", s_Molecule_MoleculeAtIndex, -1);
- rb_define_singleton_method(rb_cMolecule, "open", s_Molecule_Open, -1);
rb_define_singleton_method(rb_cMolecule, "list", s_Molecule_List, 0);
rb_define_singleton_method(rb_cMolecule, "ordered_list", s_Molecule_OrderedList, 0);
- rb_define_singleton_method(rb_cMolecule, "error_message", s_Molecule_ErrorMessage, 0);
- rb_define_singleton_method(rb_cMolecule, "set_error_message", s_Molecule_SetErrorMessage, 1);
- rb_define_singleton_method(rb_cMolecule, "error_message=", s_Molecule_SetErrorMessage, 1);
+
+ rb_define_method(rb_cMolecule, "call_subprocess_async", s_Molecule_CallSubProcessAsync, -1);
/* class MolEnumerable */
rb_cMolEnumerable = rb_define_class_under(rb_mMolby, "MolEnumerable", rb_cObject);
OSDN Git Service
