s_FractRSym, s_FractXSym, s_FractYSym, s_FractZSym,
s_SigmaSym, s_SigmaXSym, s_SigmaYSym, s_SigmaZSym,
s_VSym, s_FSym, s_OccupancySym, s_TempFactorSym,
- s_AnisoSym, s_SymopSym, s_IntChargeSym, s_FixForceSym,
- s_FixPosSym, s_ExclusionSym, s_MMExcludeSym, s_PeriodicExcludeSym,
- s_HiddenSym, s_AnchorListSym, s_UFFTypeSym;
+ s_AnisoSym, s_AnisoEigvalSym, s_SymopSym, s_IntChargeSym,
+ s_FixForceSym, s_FixPosSym, s_ExclusionSym, s_MMExcludeSym,
+ s_PeriodicExcludeSym, s_HiddenSym, s_AnchorListSym, s_UFFTypeSym;
/* Symbols for parameter attributes */
static VALUE
s_CutoffSym, s_RadiusSym, s_ColorSym, s_FullNameSym, s_VdwRadiusSym,
s_CommentSym, s_SourceSym;
+/* Symbols for graphics */
+static VALUE
+ s_LineSym, s_PolySym, s_CylinderSym, s_ConeSym, s_EllipsoidSym;
+
/*
* Utility function
* Get ary[i] by calling "[]" method
#if __WXMSW__
char *p = strdup(StringValuePtr(*valp));
translate_char(p, '/', '\\');
- *valp = rb_str_new2(p);
+ *valp = Ruby_NewEncodedStringValue2(p);
free(p);
return StringValuePtr(*valp);
#else
VALUE retval;
char *p = strdup(fstr);
translate_char(p, '\\', '/');
- retval = rb_enc_str_new(p, strlen(p), rb_default_external_encoding());
+ retval = Ruby_NewEncodedStringValue2(p);
free(p);
return retval;
#else
- return rb_str_new2(fstr);
+ return Ruby_NewEncodedStringValue2(fstr);
#endif
}
}
VALUE
+Ruby_NewEncodedStringValue2(const char *str)
+{
+ return Ruby_NewEncodedStringValue(str, -1);
+}
+
+VALUE
Ruby_ObjToStringObj(VALUE val)
{
switch (TYPE(val)) {
s_Kernel_MessageBox(int argc, VALUE *argv, VALUE self)
{
char *str, *title, *s;
- int buttons, icon;
+ int buttons, icon, retval;
VALUE sval, tval, bval, ival;
rb_scan_args(argc, argv, "22", &sval, &tval, &bval, &ival);
str = StringValuePtr(sval);
else
rb_raise(rb_eMolbyError, "the icon specification should be either :info, :warning or :error");
} else icon = 1;
- MyAppCallback_messageBox(str, title, buttons, icon);
- return Qnil;
+ retval = MyAppCallback_messageBox(str, title, buttons, icon);
+ return (retval ? Qtrue : Qfalse);
}
/*
} else buf[0] = 0;
retval = MyAppCallback_getTextWithPrompt(StringValuePtr(prompt), buf, sizeof buf);
if (retval)
- return rb_str_new2(buf);
+ return Ruby_NewEncodedStringValue2(buf);
else
return Qnil;
}
s_Kernel_ExportToClipboard(VALUE self, VALUE sval)
{
#if !defined(__CMDMAC__)
- const char *s = StringValuePtr(sval);
+ const char *s;
char *ns;
+ if (!gUseGUI)
+ return Qnil;
+ s = StringValuePtr(sval);
#if __WXMSW__
/* Convert the end-of-line characters */
{ const char *p; int nc; char *np;
/*
* call-seq:
+ * hartree_to_kcal(val)
+ *
+ * Convert hartree to kcal/mol
+ */
+static VALUE
+s_Kernel_HartreeToKcal(VALUE self, VALUE fval)
+{
+ double d = NUM2DBL(rb_Float(fval));
+ return rb_float_new(d * 627.5094740630557);
+}
+
+/*
+ * call-seq:
+ * kcal_to_hartree(val)
+ *
+ * Convert kcal/mol to hartree
+ */
+static VALUE
+s_Kernel_KcalToHartree(VALUE self, VALUE fval)
+{
+ double d = NUM2DBL(rb_Float(fval));
+ return rb_float_new(d / 627.5094740630557);
+}
+
+/*
+ * call-seq:
+ * hartree_to_kj(val)
+ *
+ * Convert hartree to kJ/mol
+ */
+static VALUE
+s_Kernel_HartreeToKJ(VALUE self, VALUE fval)
+{
+ double d = NUM2DBL(rb_Float(fval));
+ return rb_float_new(d * 2625.4996394798253);
+}
+
+/*
+ * call-seq:
+ * kj_to_hartree(val)
+ *
+ * Convert kJ/mol to hartree
+ */
+static VALUE
+s_Kernel_KJToHartree(VALUE self, VALUE fval)
+{
+ double d = NUM2DBL(rb_Float(fval));
+ return rb_float_new(d / 2625.4996394798253);
+}
+
+/*
+ * call-seq:
+ * bohr_to_angstrom(val)
+ *
+ * Convert bohr to angstrom
+ */
+static VALUE
+s_Kernel_BohrToAngstrom(VALUE self, VALUE fval)
+{
+ double d = NUM2DBL(rb_Float(fval));
+ return rb_float_new(d * 0.529177210903);
+}
+
+/*
+ * call-seq:
+ * angstrom_to_bohr(val)
+ *
+ * Convert angstrom to bohr
+ */
+static VALUE
+s_Kernel_AngstromToBohr(VALUE self, VALUE fval)
+{
+ double d = NUM2DBL(rb_Float(fval));
+ return rb_float_new(d / 0.529177210903);
+}
+
+/*
+ * call-seq:
* stdout.write(str)
*
* Put the message in the main text view in black color.
s_StandardInputGets(int argc, VALUE *argv, VALUE self)
{
VALUE pval, rval;
- pval = rb_str_new2("Enter a line:");
+ pval = Ruby_NewEncodedStringValue2("Enter a line:");
rval = s_Kernel_Ask(1, &pval, self);
if (rval == Qnil)
rb_interrupt();
oldval = s_interrupt_flag;
if (val != Qundef) {
s_interrupt_flag = val;
- if (val == Qfalse) {
- s_HideProgressPanel(self);
- }
}
return oldval;
}
sNonEmptySym = ID2SYM(rb_intern("non_empty"));
sSelectionSym = ID2SYM(rb_intern("selection"));
sMolProc = rb_eval_string("lambda { |m| m != nil }");
+ rb_define_variable("$is_a_molecule_proc", &sMolProc);
sNonEmptyProc = rb_eval_string("lambda { |m| m.is_a?(Molecule) && m.natoms > 0 }");
+ rb_define_variable("$is_molecule_not_empty_proc", &sNonEmptyProc);
sSelectionProc = rb_eval_string("lambda { |m| m.is_a?(Molecule) && m.selection.count > 0 }");
+ rb_define_variable("$has_molecule_selection_proc", &sSelectionProc);
sTrueProc = rb_eval_string("lambda { |m| true }");
- rb_global_variable(&sMolProc);
- rb_global_variable(&sNonEmptyProc);
- rb_global_variable(&sSelectionProc);
- rb_global_variable(&sTrueProc);
+ rb_define_variable("$always_true_proc", &sTrueProc);
}
if (pval == Qnil) {
s_Ruby_UpdateUI_handler(VALUE data)
{
void **p = (void **)data;
- int index = (int)p[0];
+ int index = (intptr_t)p[0];
Molecule *mol = (Molecule *)p[1];
int *outChecked = (int *)p[2];
char **outTitle = (char **)p[3];
int status;
void *p[4];
VALUE retval;
- p[0] = (void *)index;
+ p[0] = (void *)(intptr_t)index;
p[1] = mol;
p[2] = outChecked;
p[3] = outTitle;
*
* Call subprocess. A progress dialog window is displayed, with a message
* "Running #{process_name}...".
+ * cmd is either a single string of an array of string. If it is a single string, then
+ * it will be given to wxExecute as a single argument. In this case, the string can be
+ * split into arguments by whitespace. If this behavior is not intended, then use an array
+ * containing a single string.
* A callback proc can be given, which is called periodically during execution. If the proc returns
* nil or false, then the execution will be interrupted.
* If stdout_file or stderr_file is a filename, then the message will be sent to the file; if the
s_Kernel_CallSubProcess(int argc, VALUE *argv, VALUE self)
{
VALUE cmd, procname, cproc, stdout_val, stderr_val;
- int n;
+ VALUE save_interruptFlag;
+ int n, exitstatus, pid;
char *sout, *serr;
+ const char *pnamestr, **cmdargv;
FILE *fpout, *fperr;
rb_scan_args(argc, argv, "23", &cmd, &procname, &cproc, &stdout_val, &stderr_val);
rb_raise(rb_eMolbyError, "Cannot open file for standard output (%s)", serr);
}
}
+
+ save_interruptFlag = s_SetInterruptFlag(self, Qnil);
+ if (procname != Qnil)
+ pnamestr = StringValuePtr(procname);
+ else pnamestr = NULL;
+ if (rb_obj_is_kind_of(cmd, rb_cString)) {
+ cmdargv = calloc(sizeof(cmdargv[0]), 3);
+ cmdargv[0] = StringValuePtr(cmd);
+ cmdargv[1] = "";
+ cmdargv[2] = NULL;
+ } else {
+ cmd = rb_ary_to_ary(cmd);
+ cmdargv = calloc(sizeof(cmdargv[0]), RARRAY_LEN(cmd) + 1);
+ for (n = 0; n < RARRAY_LEN(cmd); n++) {
+ cmdargv[n] = StringValuePtr(RARRAY_PTR(cmd)[n]);
+ }
+ cmdargv[n] = NULL;
+ }
+ n = MyAppCallback_callSubProcess(cmdargv, pnamestr, (cproc == Qnil ? NULL : s_Kernel_CallSubProcess_Callback), (cproc == Qnil ? NULL : (void *)cproc), fpout, fperr, &exitstatus, &pid);
+ s_SetInterruptFlag(self, save_interruptFlag);
+ free(cmdargv);
- n = MyAppCallback_callSubProcess(StringValuePtr(cmd), StringValuePtr(procname), (cproc == Qnil ? NULL : s_Kernel_CallSubProcess_Callback), (cproc == Qnil ? NULL : (void *)cproc), fpout, fperr);
-
if (fpout != NULL && fpout != (FILE *)1)
fclose(fpout);
if (fperr != NULL && fperr != (FILE *)1)
s_Kernel_Backquote(VALUE self, VALUE cmd)
{
char *buf;
- int n;
+ int n, exitstatus, pid;
VALUE val;
- n = MyAppCallback_callSubProcess(StringValuePtr(cmd), NULL, DUMMY_CALLBACK, &buf, NULL, NULL);
- if (n != 0)
- rb_raise(rb_eMolbyError, "Cannot invoke command '%s'", StringValuePtr(cmd));
- if (buf != NULL) {
+ const char *cmdargv[3];
+ cmdargv[0] = StringValuePtr(cmd);
+ cmdargv[1] = "";
+ cmdargv[2] = NULL;
+ n = MyAppCallback_callSubProcess(cmdargv, NULL, DUMMY_CALLBACK, &buf, NULL, NULL, &exitstatus, &pid);
+/* fprintf(stderr, "n = %d, exitstatus = %d, pid = %d\n", n, exitstatus, pid); */
+ if (n >= 0 && buf != NULL) {
val = Ruby_NewEncodedStringValue(buf, 0);
free(buf);
} else {
val = Ruby_NewEncodedStringValue("", 0);
}
+ rb_last_status_set(exitstatus, pid);
return val;
}
#pragma mark ====== IO extension ======
+static VALUE
+s_Ruby_str_encode_protected(VALUE val)
+{
+ return rb_str_encode(val, rb_enc_from_encoding(rb_default_external_encoding()),
+ ECONV_INVALID_REPLACE | ECONV_UNDEF_REPLACE, Qnil);
+}
+
/*
* call-seq:
* gets_any_eol
static VALUE
s_IO_gets_any_eol(VALUE self)
{
- VALUE val, cval;
+ VALUE val, val2, cval;
char buf[1024];
- int i, c;
+ int i, c, status;
static ID id_getbyte = 0, id_ungetbyte;
if (id_getbyte == 0) {
id_getbyte = rb_intern("getbyte");
if (c != 0x0a)
rb_funcall(self, id_ungetbyte, 1, cval);
}
+ break;
} else if (c != 0x0a) {
buf[i++] = c;
if (i >= 1020) {
val = rb_str_new(buf, i);
else if (i > 0)
rb_str_append(val, rb_str_new(buf, i));
- val = rb_str_encode(val, rb_enc_from_encoding(rb_default_external_encoding()), 0, Qnil);
+ val2 = rb_protect(s_Ruby_str_encode_protected, val, &status); /* Ignore exception */
+ if (status == 0)
+ val = val2;
if (cval != Qnil) {
/* Needs a end-of-line mark */
cval = rb_gv_get("$/");
Int tp;
s_UnionParFromValue(self, &tp, 0);
if (tp == kElementParType)
- return rb_str_new2("element");
+ return Ruby_NewEncodedStringValue2("element");
tp -= kFirstParType;
if (tp >= 0 && tp < sizeof(s_ParameterTypeNames) / sizeof(s_ParameterTypeNames[0]))
- return rb_str_new2(s_ParameterTypeNames[tp]);
+ return Ruby_NewEncodedStringValue2(s_ParameterTypeNames[tp]);
else rb_raise(rb_eMolbyError, "Internal error: parameter type tag is out of range (%d)", tp);
}
char name[5];
strncpy(name, up->atom.name, 4);
name[4] = 0;
- return rb_str_new2(name);
+ return Ruby_NewEncodedStringValue2(name);
} else rb_raise(rb_eMolbyError, "invalid member name");
}
char fullname[16];
strncpy(fullname, up->atom.fullname, 15);
fullname[15] = 0;
- return rb_str_new2(fullname);
+ return Ruby_NewEncodedStringValue2(fullname);
} else rb_raise(rb_eMolbyError, "invalid member fullname");
}
com = up->bond.com;
if (com == 0)
return Qnil;
- else return rb_str_new2(ParameterGetComment(com));
+ else return Ruby_NewEncodedStringValue2(ParameterGetComment(com));
}
/*
return Qfalse; /* undefined */
else if (src == 0)
return Qnil; /* local */
- else return rb_str_new2(ParameterGetComment(src));
+ else return Ruby_NewEncodedStringValue2(ParameterGetComment(src));
}
static void
snprintf(buf, sizeof buf, "element %2.2s %3d %6.3f %6.3f %6.3f %6.3f %8.4f %s %6.3f", up->atom.name, up->atom.number, up->atom.radius, up->atom.red / 65535.0, up->atom.green / 65535.0, up->atom.blue / 65535.0, up->atom.weight, up->atom.fullname, up->atom.vdw_radius);
break;
}
- return rb_str_new2(buf);
+ return Ruby_NewEncodedStringValue2(buf);
}
/*
Data_Get_Struct(self, ParEnumerable, pen);
tp = pen->parType;
if (tp == kElementParType)
- return rb_str_new2("element");
+ return Ruby_NewEncodedStringValue2("element");
tp -= kFirstParType;
if (tp >= 0 && tp < sizeof(s_ParameterTypeNames) / sizeof(s_ParameterTypeNames[0]))
- return rb_str_new2(s_ParameterTypeNames[tp]);
+ return Ruby_NewEncodedStringValue2(s_ParameterTypeNames[tp]);
else rb_raise(rb_eMolbyError, "Internal error: parameter type tag is out of range (%d)", tp);
}
static VALUE s_AtomRef_GetSegName(VALUE self) {
char *p = s_AtomFromValue(self)->segName;
- return rb_str_new(p, strlen_limit(p, 4));
+ return Ruby_NewEncodedStringValue(p, strlen_limit(p, 4));
}
static VALUE s_AtomRef_GetResSeq(VALUE self) {
static VALUE s_AtomRef_GetResName(VALUE self) {
char *p = s_AtomFromValue(self)->resName;
- return rb_str_new(p, strlen_limit(p, 4));
+ return Ruby_NewEncodedStringValue(p, strlen_limit(p, 4));
}
static VALUE s_AtomRef_GetName(VALUE self) {
char *p = s_AtomFromValue(self)->aname;
- return rb_str_new(p, strlen_limit(p, 4));
+ return Ruby_NewEncodedStringValue(p, strlen_limit(p, 4));
}
static VALUE s_AtomRef_GetAtomType(VALUE self) {
int type = s_AtomFromValue(self)->type;
char *p = (type == 0 ? "" : AtomTypeDecodeToString(type, NULL));
- return rb_str_new(p, strlen_limit(p, 6));
+ return Ruby_NewEncodedStringValue(p, strlen_limit(p, 6));
}
static VALUE s_AtomRef_GetCharge(VALUE self) {
static VALUE s_AtomRef_GetElement(VALUE self) {
char *p = s_AtomFromValue(self)->element;
- return rb_str_new(p, strlen_limit(p, 4));
+ return Ruby_NewEncodedStringValue(p, strlen_limit(p, 4));
}
static VALUE s_AtomRef_GetAtomicNumber(VALUE self) {
return retval;
}
+static VALUE s_AtomRef_GetAnisoEigenValues(VALUE self) {
+ VALUE retval;
+ int i;
+ Atom *ap = s_AtomFromValue(self);
+ if (ap->aniso == NULL)
+ return Qnil;
+ retval = rb_ary_new();
+ for (i = 0; i < 3; i++)
+ rb_ary_push(retval, rb_float_new(ap->aniso->eigval[i]));
+ return retval;
+}
+
static VALUE s_AtomRef_GetSymop(VALUE self) {
VALUE retval;
Atom *ap = s_AtomFromValue(self);
static VALUE s_AtomRef_GetUFFType(VALUE self) {
char *p = s_AtomFromValue(self)->uff_type;
- return rb_str_new(p, strlen_limit(p, 5));
+ return Ruby_NewEncodedStringValue(p, strlen_limit(p, 5));
}
static VALUE s_AtomRef_SetIndex(VALUE self, VALUE val) {
return val;
}
+static VALUE s_AtomRef_SetAnisoEigenValues(VALUE self, VALUE val) {
+ rb_raise(rb_eMolbyError, "Eigenvalues of anisotropic factors are read-only.");
+ return val; /* Not reached */
+}
+
static VALUE s_AtomRef_SetSymop(VALUE self, VALUE val) {
Molecule *mol;
Atom *ap;
{"occupancy", &s_OccupancySym, 0, s_AtomRef_GetOccupancy, s_AtomRef_SetOccupancy},
{"temp_factor", &s_TempFactorSym, 0, s_AtomRef_GetTempFactor, s_AtomRef_SetTempFactor},
{"aniso", &s_AnisoSym, 0, s_AtomRef_GetAniso, s_AtomRef_SetAniso},
+ {"aniso_eigenvalues", &s_AnisoEigvalSym, 0, s_AtomRef_GetAnisoEigenValues, s_AtomRef_SetAnisoEigenValues},
{"symop", &s_SymopSym, 0, s_AtomRef_GetSymop, s_AtomRef_SetSymop},
{"int_charge", &s_IntChargeSym, 0, s_AtomRef_GetIntCharge, s_AtomRef_SetIntCharge},
{"fix_force", &s_FixForceSym, 0, s_AtomRef_GetFixForce, s_AtomRef_SetFixForce},
if (idx2 < 0 || idx2 >= mol->nresidues)
rb_raise(rb_eIndexError, "residue index out of range (%d; should be %d..%d)", idx1, -mol->nresidues, mol->nresidues - 1);
p = mol->residues[idx2];
- return rb_str_new(p, strlen_limit(p, 4));
+ return Ruby_NewEncodedStringValue(p, strlen_limit(p, 4));
}
}
return Qnil;
#pragma mark ====== Molecule Class ======
-/* An malloc'ed string buffer. Retains the error/warning message from the last ***load/***save method. */
+#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;
s_Molecule_AtomGroupFromValue(VALUE self, VALUE val)
{
IntGroup *ig;
+ Molecule *mp1;
+ Data_Get_Struct(self, Molecule, mp1);
val = rb_funcall(self, rb_intern("atom_group"), 1, val);
if (!rb_obj_is_kind_of(val, rb_cIntGroup))
rb_raise(rb_eMolbyError, "IntGroup instance is expected");
Data_Get_Struct(val, IntGroup, ig);
+ IntGroupRemove(ig, mp1->natoms, ATOMS_MAX_NUMBER); /* Limit the group member to existing atoms */
IntGroupRetain(ig);
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, int isloading, const char *msg, const char *fname)
+{
+ if (gLoadSaveErrorMessage != NULL) {
+ MyAppCallback_setConsoleColor(1);
+ MyAppCallback_showScriptMessage("On %s %s:\n%s\n", (isloading ? "loading" : "saving"), 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.
rb_scan_args(argc, argv, "1", &fname);
fstr = FileStringValuePtr(fname);
retval = MoleculeLoadMbsfFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load mbsf", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load mbsf", fstr);
return Qtrue;
}
rb_scan_args(argc, argv, "11", &fname, &pdbname);
fstr = FileStringValuePtr(fname);
retval = MoleculeLoadPsfFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load psf", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load psf", fstr);
pdbstr = NULL;
if (!NIL_P(pdbname)) {
pdbstr = strdup(FileStringValuePtr(pdbname));
retval = MoleculeReadCoordinatesFromPdbFile(mol, pdbstr, &gLoadSaveErrorMessage);
free(pdbstr);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load coordinates from pdb", pdbstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load coordinates from pdb", pdbstr);
}
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeReadCoordinatesFromPdbFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load pdb", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load pdb", fstr);
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeReadCoordinatesFromDcdFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load dcd", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load dcd", fstr);
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeLoadTepFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load ORTEP file", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load ORTEP file", fstr);
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeLoadShelxFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load SHELX res file", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load SHELX res file", fstr);
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeLoadGaussianFchkFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load Gaussian fchk", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load Gaussian fchk", fstr);
return Qtrue;
}
MyAppCallback_showProgressPanel("Loading GAMESS dat file...");
retval = MoleculeLoadGamessDatFile(mol, fstr, &gLoadSaveErrorMessage);
MyAppCallback_hideProgressPanel();
- s_Molecule_RaiseOnLoadSave(retval, "Failed to load GAMESS dat", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 1, "Failed to load GAMESS dat", fstr);
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeWriteToMbsfFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to save mbsf", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 0, "Failed to save mbsf", fstr);
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeWriteToPsfFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to save psf", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 0, "Failed to save psf", fstr);
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeWriteToPdbFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to save pdb", fstr);
+ s_Molecule_RaiseOnLoadSave(retval, 0, "Failed to save pdb", fstr);
return Qtrue;
}
MoleculeClearLoadSaveErrorMessage();
fstr = FileStringValuePtr(fname);
retval = MoleculeWriteToDcdFile(mol, fstr, &gLoadSaveErrorMessage);
- s_Molecule_RaiseOnLoadSave(retval, "Failed to save dcd", fstr);
- 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);
+ s_Molecule_RaiseOnLoadSave(retval, 0, "Failed to save dcd", fstr);
return Qtrue;
}
-*/
/* load([ftype, ] fname, ...) */
static VALUE
failure:
rval = rb_str_to_str(argv[0]);
asprintf(&p, "Failed to %s file %s", (loadFlag ? "load" : "save"), type);
- s_Molecule_RaiseOnLoadSave(1, p, StringValuePtr(rval));
+ s_Molecule_RaiseOnLoadSave(1, loadFlag, p, StringValuePtr(rval));
return Qnil; /* Does not reach here */
success:
/*
* 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;
+
+ if (!gUseGUI) {
+ rb_raise(rb_eMolbyError, "Molecule.open is not usable in non-GUI mode. Use Molecule.new instead.");
+ }
+
+ 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 Ruby_NewEncodedStringValue2(gLoadSaveErrorMessage);
+}
+
+/*
+ * call-seq:
+ * set_error_message(String)
+ * Molecule.error_message = String
+ *
+ * Set the error_message for the present load/save method.
+ */
+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:
+ * set_molecule(Molecule)
+ *
+ * Duplicate the given molecule and set to self. The present molecule must be empty.
+ * This method is exclusively used for associating a new document with an existing molecule.
+ */
+static VALUE
+s_Molecule_SetMolecule(VALUE self, VALUE mval)
+{
+ Molecule *mp1, *mp2;
+ Data_Get_Struct(self, Molecule, mp1);
+ mp2 = MoleculeFromValue(mval);
+ MoleculeInitWithMolecule(mp1, mp2);
+ MoleculeCallback_notifyModification(mp1, 1);
+ return self;
+}
+
+#pragma mark ------ Name attributes ------
+
+/*
+ * call-seq:
* name -> String
*
* Returns the display name of the molecule. If the molecule has no associated
if (buf[0] == 0)
return Qnil;
else
- return rb_str_new2(buf);
+ return Ruby_NewEncodedStringValue2(buf);
}
/*
p = strrchr(buf, '/');
if (p != NULL)
*p = 0;
- return rb_str_new2(buf);
+ return Ruby_NewEncodedStringValue2(buf);
}
}
if (buf[0] == 0) {
/* No associated document */
snprintf(buf, sizeof buf, "#<Molecule:0x%lx>", self);
- return rb_str_new2(buf);
+ return Ruby_NewEncodedStringValue2(buf);
} else {
/* Check whether the document name is duplicate */
char buf2[256];
} else {
snprintf(buf2, sizeof buf2, "Molecule[\"%s\"]", buf);
}
- return rb_str_new2(buf2);
+ return Ruby_NewEncodedStringValue2(buf2);
}
}
+#pragma mark ------ MolEnumerables ------
+
+static VALUE
+s_Molecule_MolEnumerable(VALUE self, int kind)
+{
+ Molecule *mol;
+ MolEnumerable *mseq;
+ Data_Get_Struct(self, Molecule, mol);
+ mseq = MolEnumerableNew(mol, kind);
+ return Data_Wrap_Struct(rb_cMolEnumerable, 0, (void (*)(void *))MolEnumerableRelease, mseq);
+}
+
/*
* call-seq:
- * open -> Molecule
- * open(file) -> Molecule
+ * atoms -> MolEnumerable
*
- * Create a new molecule from file as a document. If file is not given, an untitled document is created.
+ * Returns a MolEnumerable object representing the array of atoms.
*/
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 */
-}
-
-static VALUE
-s_Molecule_MolEnumerable(VALUE self, int kind)
-{
- Molecule *mol;
- MolEnumerable *mseq;
- Data_Get_Struct(self, Molecule, mol);
- mseq = MolEnumerableNew(mol, kind);
- return Data_Wrap_Struct(rb_cMolEnumerable, 0, (void (*)(void *))MolEnumerableRelease, mseq);
-}
-
-/*
- * call-seq:
- * atoms -> MolEnumerable
- *
- * Returns a MolEnumerable object representing the array of atoms.
- */
-static VALUE
-s_Molecule_Atoms(VALUE self)
+s_Molecule_Atoms(VALUE self)
{
return s_Molecule_MolEnumerable(self, kAtomKind);
}
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
- *
- * Returns the MD parameter for the idx-th bond.
- */
-/*
-static VALUE
-s_Molecule_BondPar(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);
-}
-*/
-
-/*
- * call-seq:
- * angle_par(idx) -> ParameterRef
- *
- * Returns the MD parameter for the idx-th angle.
- */
-/*
-static VALUE
-s_Molecule_AnglePar(VALUE self, VALUE val)
-{
- 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);
-}
-*/
-/*
- * call-seq:
- * dihedral_par(idx) -> ParameterRef
- *
- * Returns the MD parameter for the idx-th dihedral.
- */
-/*
-static VALUE
-s_Molecule_DihedralPar(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->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);
-}
-*/
-/*
- * call-seq:
- * improper_par(idx) -> ParameterRef
- *
- * 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).
- */
-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;
-}
-
-/*
- * call-seq:
- * find_angles -> Integer
- *
- * Find the angles from the bonds. Returns the number of angles newly created.
- */
-/*
-static VALUE
-s_Molecule_FindAngles(VALUE self)
-{
- 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);
-}
-*/
-/*
- * call-seq:
- * find_dihedrals -> Integer
- *
- * Find the dihedrals from the bonds. Returns the number of dihedrals newly created.
- */
-/*
-static VALUE
-s_Molecule_FindDihedrals(VALUE self)
-{
- 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);
- }
- return INT2NUM(nip);
-}
-*/
-
/*
* call-seq:
* nresidues = Integer
return self;
}
+#pragma mark ------ Atom Group ------
+
+static VALUE
+s_Molecule_AtomGroup_i(VALUE arg, VALUE values)
+{
+ 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:
- * cell -> [a, b, c, alpha, beta, gamma [, sig_a, sig_b, sig_c, sig_alpha, sig_beta, sig_gamma]]
+ * atom_group
+ * atom_group {|aref| ...}
+ * atom_group(arg1, arg2, ...)
+ * atom_group(arg1, arg2, ...) {|aref| ...}
*
- * Returns the unit cell parameters. If cell is not set, returns nil.
- */
-static VALUE
-s_Molecule_Cell(VALUE self)
-{
- Molecule *mol;
- int i;
- VALUE val;
+ * 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.
+ *
+ */
+static VALUE
+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);
- 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]));
+ 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++;
}
}
- return val;
+ 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:
- * 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)
+ * selection -> IntGroup
*
- * 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.
+ * Returns the current selection.
*/
static VALUE
-s_Molecule_SetCell(int argc, VALUE *argv, VALUE self)
+s_Molecule_Selection(VALUE self)
{
Molecule *mol;
- VALUE val, cval;
- int i, convert_coord, n;
- double d[12];
+ IntGroup *ig;
+ VALUE val;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "11", &val, &cval);
- if (val == Qnil) {
- n = 0;
+ 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 {
- 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]));
+ val = IntGroup_Alloc(rb_cIntGroup);
}
- convert_coord = (RTEST(cval) ? 1 : 0);
- MolActionCreateAndPerform(mol, gMolActionSetCell, n, d, convert_coord);
return val;
}
-/*
- * call-seq:
- * box -> [avec, bvec, cvec, origin, flags]
- *
- * 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_Box(VALUE self)
+s_Molecule_SetSelectionSub(VALUE self, VALUE val, int undoable)
{
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);
+ 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:
- * 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
+ * selection = IntGroup
*
- * 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.
+ * 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_SetBox(VALUE self, VALUE aval)
+s_Molecule_SetSelection(VALUE self, VALUE val)
{
- 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;
- 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;
- } 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;
+ return s_Molecule_SetSelectionSub(self, val, 0);
}
/*
* call-seq:
- * cell_periodicity -> [n1, n2, n3]
+ * set_undoable_selection(IntGroup)
*
- * Get flags denoting whether the cell is periodic along the a/b/c axes. If the cell is not defined
- * nil is returned.
+ * Set the current selection with undo registration. The right-hand operand may be nil.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_CellPeriodicity(VALUE self)
+s_Molecule_SetUndoableSelection(VALUE self, VALUE val)
{
- Molecule *mol;
- 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]));
+ return s_Molecule_SetSelectionSub(self, val, 1);
}
+#pragma mark ------ Editing ------
+
/*
* call-seq:
- * self.cell_periodicity = [n1, n2, n3] or Integer or nil
- * set_cell_periodicity = [n1, n2, n3] or Integer or nil
+ * extract(group, dummy_flag = nil) -> 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.
- * This operation is undoable.
+ * 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_SetCellPeriodicity(VALUE self, VALUE arg)
+s_Molecule_Extract(int argc, VALUE *argv, VALUE self)
{
- 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));
- }
+ 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);
}
- MolActionCreateAndPerform(mol, gMolActionSetCellPeriodicity, flag);
- return arg;
+ IntGroupRelease(ig);
+ return retval;
}
/*
* call-seq:
- * cell_flexibility -> bool
+ * add(molecule2) -> self
*
- * Returns the unit cell is flexible or not
+ * 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_CellFlexibility(VALUE self)
+s_Molecule_Add(VALUE self, VALUE val)
{
- 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; */
+ Molecule *mol1, *mol2;
+ Data_Get_Struct(self, Molecule, mol1);
+ mol2 = MoleculeFromValue(val);
+ MolActionCreateAndPerform(mol1, gMolActionMergeMolecule, mol2, NULL);
+ return self;
}
/*
* call-seq:
- * self.cell_flexibility = bool
- * set_cell_flexibility(bool)
+ * remove(group) -> Molecule
*
- * Change the unit cell is flexible or not
+ * The atoms designated by the given group are removed from the molecule.
+ * This operation is undoable.
*/
static VALUE
-s_Molecule_SetCellFlexibility(VALUE self, VALUE arg)
+s_Molecule_Remove(VALUE self, VALUE group)
{
- rb_warn("set_cell_flexibility is obsolete (unit cell is always frame dependent)");
+ Molecule *mol1;
+ IntGroup *ig, *bg;
+ Int i;
+ IntGroupIterator iter;
+
+ ig = s_Molecule_AtomGroupFromValue(self, group);
+/* 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); */
+ Data_Get_Struct(self, Molecule, mol1);
+
+ /* 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;
-/* Molecule *mol;
- Data_Get_Struct(self, Molecule, mol);
- MolActionCreateAndPerform(mol, gMolActionSetCellFlexibility, RTEST(arg) != 0);
- return self; */
}
/*
* call-seq:
- * cell_transform -> Transform
+ * create_atom(name, pos = -1) -> AtomRef
*
- * Get the transform matrix that converts internal coordinates to cartesian coordinates.
- * If cell is not defined, 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_CellTransform(VALUE self)
+s_Molecule_CreateAnAtom(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;
Data_Get_Struct(self, Molecule, mol);
- if (mol == NULL || mol->cell == NULL)
+ 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;
- return ValueFromTransform(&(mol->cell->tr));
+ aref = AtomRefNew(mol, pos);
+ return Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
}
/*
* call-seq:
- * symmetry -> Array of Transforms
- * symmetries -> Array of Transforms
+ * duplicate_atom(atomref, pos = -1) -> AtomRef
*
- * Get the currently defined symmetry operations. If no symmetry operation is defined,
- * returns an empty array.
+ * 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_Symmetry(VALUE self)
+s_Molecule_DuplicateAnAtom(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE val;
- int i;
+ const Atom *apsrc;
+ Atom arec;
+ AtomRef *aref;
+ VALUE retval, aval, ival;
+ Int pos;
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]));
+ 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);
}
- return val;
+ 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:
- * nsymmetries -> Integer
+ * create_bond(n1, n2, ...) -> Integer
*
- * Get the number of currently defined symmetry operations.
+ * 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_Nsymmetries(VALUE self)
+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);
- return INT2NUM(mol->nsyms);
+ 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:
- * add_symmetry(Transform) -> Integer
+ * molecule.remove_bonds(n1, n2, ...) -> Integer
*
- * 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.
+ * 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_AddSymmetry(VALUE self, VALUE trans)
+s_Molecule_RemoveBond(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- Transform tr;
+ 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);
- TransformFromValue(trans, &tr);
- MolActionCreateAndPerform(mol, gMolActionAddSymmetryOperation, &tr);
- return INT2NUM(mol->nsyms);
-}
-
-/*
- * call-seq:
- * remove_symmetry(count = nil) -> Integer
- * remove_symmetries(count = nil) -> Integer
- *
- * 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_RemoveSymmetry(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- VALUE cval;
- int i, n;
- 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;
+ 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);
+ }
+ }
}
- for (i = 0; i < n; i++)
- MolActionCreateAndPerform(mol, gMolActionDeleteSymmetryOperation);
- return INT2NUM(mol->nsyms);
-}
-
-static VALUE
-s_Molecule_AtomGroup_i(VALUE arg, VALUE values)
-{
- 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;
+ if (bg != NULL) {
+ MolActionCreateAndPerform(mol, gMolActionDeleteBonds, bg);
+ i = IntGroupGetCount(bg);
+ IntGroupRelease(bg);
+ } else i = 0;
+ return INT2NUM(i);
}
/*
* 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.
+ * 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_AtomGroup(int argc, VALUE *argv, VALUE self)
+s_Molecule_AssignBondOrder(VALUE self, VALUE idxval, VALUE dval)
{
- IntGroup *ig1, *ig2;
Molecule *mol;
- Int i, startPt, interval;
- VALUE retval = IntGroup_Alloc(rb_cIntGroup);
- Data_Get_Struct(retval, IntGroup, ig1);
+ IntGroup *ig;
Data_Get_Struct(self, Molecule, mol);
- if (argc == 0) {
- IntGroup_RaiseIfError(IntGroupAdd(ig1, 0, mol->natoms));
+ 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 {
- 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);
+ 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]));
}
- IntGroupRelease(ig2);
+ MolActionCreateAndPerform(mol, gMolActionAssignBondOrders, n, dp, ig);
+ free(dp);
+ IntGroupRelease(ig);
}
-
- /* Remove points that are out of bounds */
- IntGroup_RaiseIfError(IntGroupRemove(ig1, mol->natoms, INT_MAX));
-
- return retval;
+ return self;
}
/*
* call-seq:
- * atom_index(val) -> Integer
+ * get_bond_order(idx) -> Float
+ * get_bond_orders(group) -> Array
*
- * 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 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_AtomIndex(VALUE self, VALUE val)
+s_Molecule_GetBondOrder(VALUE self, VALUE idxval)
{
Molecule *mol;
+ IntGroup *ig;
+ Double *dp;
+ VALUE retval;
+ Int i, n, numericArg;
Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(s_Molecule_AtomIndexFromValue(mol, val));
+ 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:
- * extract(group, dummy_flag = nil) -> Molecule
+ * bond_exist?(idx1, idx2) -> bool
*
- * 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.
+ * 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_Extract(int argc, VALUE *argv, VALUE self)
+s_Molecule_BondExist(VALUE self, VALUE ival1, VALUE ival2)
{
- 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);
+ Molecule *mol;
+ Int idx1, idx2, i;
+ 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;
}
- IntGroupRelease(ig);
- return retval;
+ return Qfalse;
}
/*
* call-seq:
- * add(molecule2) -> self
+ * add_angle(n1, n2, n3) -> Molecule
*
- * Combine two molecules. The residue numbers of the newly added atoms may be renumbered to avoid
- conflicts.
- This operation is undoable.
+ * 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_Add(VALUE self, VALUE val)
+s_Molecule_AddAngle(VALUE self, VALUE v1, VALUE v2, VALUE v3)
{
- Molecule *mol1, *mol2;
- Data_Get_Struct(self, Molecule, mol1);
- mol2 = MoleculeFromValue(val);
- MolActionCreateAndPerform(mol1, gMolActionMergeMolecule, mol2, NULL);
- return self;
+ 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(group) -> Molecule
+ * remove_angle(n1, n2, n3) -> Molecule
*
- * The atoms designated by the given group are removed from the 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_Remove(VALUE self, VALUE group)
+s_Molecule_RemoveAngle(VALUE self, VALUE v1, VALUE v2, VALUE v3)
{
- 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;
+ 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;
}
/*
* call-seq:
- * create_atom(name, pos = -1) -> AtomRef
+ * add_dihedral(n1, n2, n3, n4) -> Molecule
*
- * 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.
+ * 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_CreateAnAtom(int argc, VALUE *argv, VALUE self)
+s_Molecule_AddDihedral(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
{
+ Int n[5];
Molecule *mol;
- 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, "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);
+ 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:
- * duplicate_atom(atomref, pos = -1) -> AtomRef
+ * remove_dihedral(n1, n2, n3, n4) -> Molecule
*
- * Create a new atom with the same attributes (but no bonding information)
- * with the specified atom. Returns the reference to the new atom.
+ * 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_DuplicateAnAtom(int argc, VALUE *argv, VALUE self)
+s_Molecule_RemoveDihedral(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
{
+ Int n[5];
Molecule *mol;
- const Atom *apsrc;
- Atom arec;
- AtomRef *aref;
- VALUE retval, aval, ival;
- Int pos;
+ IntGroup *ig;
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;
+ 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:
- * create_bond(n1, n2, ...) -> Integer
+ * add_improper(n1, n2, n3, n4) -> Molecule
*
- * 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.
+ * 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_CreateBond(int argc, VALUE *argv, VALUE self)
+s_Molecule_AddImproper(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
{
+ Int n[5];
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);
+ 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:
- * molecule.remove_bonds(n1, n2, ...) -> Integer
+ * remove_improper(n1, n2, n3, n4) -> Molecule
+ * remove_improper(intgroup) -> Molecule
*
- * 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.
+ * 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_RemoveBond(int argc, VALUE *argv, VALUE self)
+s_Molecule_RemoveImproper(int argc, VALUE *argv, VALUE self)
{
+ Int n[5];
+ VALUE v1, v2, v3, v4;
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");
+ IntGroup *ig;
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 (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);
}
- if (bg != NULL) {
- MolActionCreateAndPerform(mol, gMolActionDeleteBonds, bg);
- i = IntGroupGetCount(bg);
- IntGroupRelease(bg);
- } else i = 0;
- return INT2NUM(i);
+ MolActionCreateAndPerform(mol, gMolActionDeleteImpropers, ig);
+ IntGroupRelease(ig);
+ return self;
}
/*
* call-seq:
- * assign_bond_order(idx, d1)
- * assign_bond_orders(group, [d1, d2, ...])
+ * assign_residue(group, res) -> Molecule
*
- * 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)
+ * 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_AssignBondOrder(VALUE self, VALUE idxval, VALUE dval)
+s_Molecule_AssignResidue(VALUE self, VALUE range, VALUE res)
{
Molecule *mol;
IntGroup *ig;
+ char *p, *pp, buf[16];
+ Int resid, n;
+ Atom *ap;
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);
+
+ /* 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 {
- 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]));
+ p = StringValuePtr(res);
+ pp = strchr(p, '.');
+ if (pp != NULL) {
+ resid = atoi(pp + 1);
+ n = pp - p;
+ } else {
+ resid = -1;
+ n = strlen(p);
}
- MolActionCreateAndPerform(mol, gMolActionAssignBondOrders, n, dp, ig);
- free(dp);
- IntGroupRelease(ig);
+ 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:
- * get_bond_order(idx) -> Float
- * get_bond_orders(group) -> Array
+ * offset_residue(group, offset) -> Molecule
*
- * 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).
+ * 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_GetBondOrder(VALUE self, VALUE idxval)
+s_Molecule_OffsetResidue(VALUE self, VALUE range, VALUE offset)
{
Molecule *mol;
IntGroup *ig;
- Double *dp;
- VALUE retval;
- Int i, n, numericArg;
+ int ofs, result;
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);
+ 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 retval;
+ return self;
}
/*
* call-seq:
- * bond_exist?(idx1, idx2) -> bool
+ * renumber_atoms(array) -> IntGroup
*
- * 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.
+ * 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_BondExist(VALUE self, VALUE ival1, VALUE ival2)
+s_Molecule_RenumberAtoms(VALUE self, VALUE array)
{
- Molecule *mol;
- Int idx1, idx2, i;
- Atom *ap;
- Int *cp;
+ Molecule *mol;
+ Int *new2old;
+ IntGroup *ig;
+ int i, n;
+ VALUE *valp, retval;
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;
+ 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:
- * add_angle(n1, n2, n3) -> Molecule
+ * set_atom_attr(index, key, value)
*
- * 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.
+ * Set the atom attribute for the specified atom.
* This operation is undoable.
*/
static VALUE
-s_Molecule_AddAngle(VALUE self, VALUE v1, VALUE v2, VALUE v3)
+s_Molecule_SetAtomAttr(VALUE self, VALUE idx, VALUE key, VALUE val)
{
- Int n[4];
- Molecule *mol;
+ Molecule *mol;
+ VALUE aref, oldval;
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;
+ 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:
- * remove_angle(n1, n2, n3) -> Molecule
+ * get_atom_attr(index, key)
*
- * 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.
+ * Get the atom attribute for the specified atom.
*/
static VALUE
-s_Molecule_RemoveAngle(VALUE self, VALUE v1, VALUE v2, VALUE v3)
+s_Molecule_GetAtomAttr(VALUE self, VALUE idx, VALUE key)
{
- 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;
+ return s_Molecule_SetAtomAttr(self, idx, key, Qundef);
}
+#pragma mark ------ Undo Support ------
+
/*
* call-seq:
- * add_dihedral(n1, n2, n3, n4) -> Molecule
+ * register_undo(script, *args)
*
- * 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.
+ * Register an undo operation with the current molecule.
*/
static VALUE
-s_Molecule_AddDihedral(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
+s_Molecule_RegisterUndo(int argc, VALUE *argv, VALUE self)
{
- Int n[5];
- Molecule *mol;
+ Molecule *mol;
+ VALUE script, args;
+ MolAction *act;
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;
+ rb_scan_args(argc, argv, "1*", &script, &args);
+ act = MolActionNew(SCRIPT_ACTION("R"), StringValuePtr(script), args);
+ MolActionCallback_registerUndo(mol, act);
+ return script;
}
/*
* call-seq:
- * remove_dihedral(n1, n2, n3, n4) -> Molecule
+ * undo_enabled? -> bool
*
- * 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.
+ * Returns true if undo is enabled for this molecule; otherwise no.
*/
static VALUE
-s_Molecule_RemoveDihedral(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
+s_Molecule_UndoEnabled(VALUE self)
{
- Int n[5];
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);
- 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;
+ if (MolActionCallback_isUndoRegistrationEnabled(mol))
+ return Qtrue;
+ else return Qfalse;
}
/*
* call-seq:
- * add_improper(n1, n2, n3, n4) -> Molecule
+ * undo_enabled = bool
*
- * 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.
+ * Enable or disable undo.
*/
static VALUE
-s_Molecule_AddImproper(VALUE self, VALUE v1, VALUE v2, VALUE v3, VALUE v4)
+s_Molecule_SetUndoEnabled(VALUE self, VALUE val)
{
- 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 (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;
+ MolActionCallback_setUndoRegistrationEnabled(mol, (val != Qfalse && val != Qnil));
+ return val;
+}
+
+#pragma mark ------ Measure ------
+
+static void
+s_Molecule_DoCenterOfMass(Molecule *mol, Vector *outv, IntGroup *ig)
+{
+ switch (MoleculeCenterOfMass(mol, outv, ig)) {
+ case 2: rb_raise(rb_eMolbyError, "atom group is empty"); break;
+ case 3: rb_raise(rb_eMolbyError, "weight is zero --- atomic weights are not defined?"); break;
+ case 0: break;
+ default: rb_raise(rb_eMolbyError, "cannot calculate center of mass"); break;
+ }
}
/*
* call-seq:
- * remove_improper(n1, n2, n3, n4) -> Molecule
- * remove_improper(intgroup) -> Molecule
+ * center_of_mass(group = nil) -> Vector3D
*
- * 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 center of mass for the given set of atoms. The argument
+ * group is null, then all atoms are considered.
*/
static VALUE
-s_Molecule_RemoveImproper(int argc, VALUE *argv, VALUE self)
+s_Molecule_CenterOfMass(int argc, VALUE *argv, VALUE self)
{
- Int n[5];
- VALUE v1, v2, v3, v4;
Molecule *mol;
+ VALUE group;
IntGroup *ig;
+ Vector v;
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);
+ rb_scan_args(argc, argv, "01", &group);
+ ig = (NIL_P(group) ? NULL : s_Molecule_AtomGroupFromValue(self, group));
+ s_Molecule_DoCenterOfMass(mol, &v, ig);
+ if (ig != NULL)
+ IntGroupRelease(ig);
+ return ValueFromVector(&v);
+}
+
+/*
+ * call-seq:
+ * centralize(group = nil) -> self
+ *
+ * Translate the molecule so that the center of mass of the given group is located
+ * at (0, 0, 0). Equivalent to molecule.translate(molecule.center_of_mass(group) * -1).
+ */
+static VALUE
+s_Molecule_Centralize(int argc, VALUE *argv, VALUE self)
+{
+ Molecule *mol;
+ VALUE group;
+ IntGroup *ig;
+ Vector v;
+ Data_Get_Struct(self, Molecule, mol);
+ rb_scan_args(argc, argv, "01", &group);
+ ig = (NIL_P(group) ? NULL : s_Molecule_AtomGroupFromValue(self, group));
+ s_Molecule_DoCenterOfMass(mol, &v, ig);
+ if (ig != NULL)
+ IntGroupRelease(ig);
+ v.x = -v.x;
+ v.y = -v.y;
+ v.z = -v.z;
+ MolActionCreateAndPerform(mol, gMolActionTranslateAtoms, &v, NULL);
return self;
}
/*
* call-seq:
- * assign_residue(group, res) -> Molecule
+ * bounds(group = nil) -> [min, max]
*
- * 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.
+ * Calculate the boundary. The return value is an array of two Vector3D objects.
*/
static VALUE
-s_Molecule_AssignResidue(VALUE self, VALUE range, VALUE res)
+s_Molecule_Bounds(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
+ VALUE group;
IntGroup *ig;
- char *p, *pp, buf[16];
- Int resid, n;
+ Vector vmin, vmax;
+ int n;
Atom *ap;
Data_Get_Struct(self, Molecule, mol);
-
- /* 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;
+ 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;
}
- ig = s_Molecule_AtomGroupFromValue(self, range);
- if (ig == NULL || IntGroupGetCount(ig) == 0)
- return Qnil;
+ return rb_ary_new3(2, ValueFromVector(&vmin), ValueFromVector(&vmax));
+}
- 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);
+/* Get atom position or a vector */
+static void
+s_Molecule_GetVectorFromArg(Molecule *mol, VALUE val, Vector *vp)
+{
+ if (rb_obj_is_kind_of(val, rb_cInteger) || rb_obj_is_kind_of(val, rb_cString)) {
+ int n1 = s_Molecule_AtomIndexFromValue(mol, val);
+ *vp = ATOM_AT_INDEX(mol->atoms, n1)->r;
+ } else {
+ VectorFromValue(val, vp);
}
- IntGroupRelease(ig);
- return self;
}
/*
* call-seq:
- * offset_residue(group, offset) -> Molecule
+ * measure_bond(n1, n2) -> Float
*
- * Offset the residue number of the specified atoms. If any of the residue number gets
- * negative, then exception is thrown.
- * 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_OffsetResidue(VALUE self, VALUE range, VALUE offset)
+s_Molecule_MeasureBond(VALUE self, VALUE nval1, VALUE nval2)
{
Molecule *mol;
- IntGroup *ig;
- int ofs, result;
+ Vector v1, v2;
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;
+ s_Molecule_GetVectorFromArg(mol, nval1, &v1);
+ s_Molecule_GetVectorFromArg(mol, nval2, &v2);
+ return rb_float_new(MoleculeMeasureBond(mol, &v1, &v2));
}
/*
* call-seq:
- * renumber_atoms(array) -> IntGroup
+ * measure_angle(n1, n2, n3) -> Float
*
- * 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.
+ * 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_RenumberAtoms(VALUE self, VALUE array)
+s_Molecule_MeasureAngle(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3)
{
Molecule *mol;
- Int *new2old;
- IntGroup *ig;
- int i, n;
- VALUE *valp, retval;
+ Vector v1, v2, v3;
+ Double d;
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;
+ s_Molecule_GetVectorFromArg(mol, nval1, &v1);
+ s_Molecule_GetVectorFromArg(mol, nval2, &v2);
+ s_Molecule_GetVectorFromArg(mol, nval3, &v3);
+ d = MoleculeMeasureAngle(mol, &v1, &v2, &v3);
+ if (isnan(d))
+ return Qnil; /* Cannot define */
+ else return rb_float_new(d);
+}
+
+/*
+ * call-seq:
+ * measure_dihedral(n1, n2, n3, n4) -> Float
+ *
+ * Calculate the dihedral angle. The arguments can either be atom indices, the "residue:name" representation,
+ * or Vector3D values. The return value is in degree.
+ * If the crystallographic cell is defined, then the internal coordinates are convereted to the cartesian.
+ */
+static VALUE
+s_Molecule_MeasureDihedral(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3, VALUE nval4)
+{
+ Molecule *mol;
+ Vector v1, v2, v3, v4;
+ Double d;
+ Data_Get_Struct(self, Molecule, mol);
+ s_Molecule_GetVectorFromArg(mol, nval1, &v1);
+ s_Molecule_GetVectorFromArg(mol, nval2, &v2);
+ s_Molecule_GetVectorFromArg(mol, nval3, &v3);
+ s_Molecule_GetVectorFromArg(mol, nval4, &v4);
+ d = MoleculeMeasureDihedral(mol, &v1, &v2, &v3, &v4);
+ if (isnan(d))
+ return Qnil; /* Cannot define */
+ else return rb_float_new(d);
+}
+
+/*
+ * call-seq:
+ * find_conflicts(limit[, group1[, group2 [, ignore_exclusion]]]) -> [[n1, n2], [n3, n4], ...]
+ *
+ * Find pairs of atoms that are within the limit distance. If group1 and group2 are given, the
+ * first and second atom in the pair should belong to group1 and group2, respectively.
+ * If ignore_exclusion is true, then 1-2 (bonded), 1-3, 1-4 pairs are also considered.
+ */
+static VALUE
+s_Molecule_FindConflicts(int argc, VALUE *argv, VALUE self)
+{
+ Molecule *mol;
+ VALUE limval, gval1, gval2, rval, igval;
+ IntGroup *ig1, *ig2;
+ IntGroupIterator iter1, iter2;
+ Int npairs, *pairs;
+ Int n[2], i;
+ Double lim;
+ Vector r1;
+ Atom *ap1, *ap2;
+ MDExclusion *exinfo;
+ Int *exlist;
+
+ Data_Get_Struct(self, Molecule, mol);
+ rb_scan_args(argc, argv, "13", &limval, &gval1, &gval2, &igval);
+ lim = NUM2DBL(rb_Float(limval));
+ if (lim <= 0.0)
+ rb_raise(rb_eMolbyError, "the limit (%g) should be positive", lim);
+ if (gval1 != Qnil)
+ ig1 = s_Molecule_AtomGroupFromValue(self, gval1);
+ else
+ ig1 = IntGroupNewWithPoints(0, mol->natoms, -1);
+ if (gval2 != Qnil)
+ ig2 = s_Molecule_AtomGroupFromValue(self, gval2);
+ else
+ ig2 = IntGroupNewWithPoints(0, mol->natoms, -1);
+
+ if (!RTEST(igval)) {
+ /* Use the exclusion table in MDArena */
+ if (mol->par == NULL || mol->arena == NULL || mol->arena->is_initialized == 0 || mol->needsMDRebuild) {
+ VALUE mval = ValueFromMolecule(mol);
+ s_RebuildMDParameterIfNecessary(mval, Qnil);
+ }
+ exinfo = mol->arena->exinfo; /* May be NULL */
+ exlist = mol->arena->exlist;
+ } else {
+ exinfo = NULL;
+ exlist = NULL;
+ }
+ IntGroupIteratorInit(ig1, &iter1);
+ IntGroupIteratorInit(ig2, &iter2);
+ npairs = 0;
+ pairs = NULL;
+ while ((n[0] = IntGroupIteratorNext(&iter1)) >= 0) {
+ Int exn1, exn2;
+ ap1 = ATOM_AT_INDEX(mol->atoms, n[0]);
+ r1 = ap1->r;
+ if (exinfo != NULL) {
+ exn1 = exinfo[n[0]].index1;
+ exn2 = exinfo[n[0] + 1].index1;
+ } else exn1 = exn2 = -1;
+ IntGroupIteratorReset(&iter2);
+ while ((n[1] = IntGroupIteratorNext(&iter2)) >= 0) {
+ ap2 = ATOM_AT_INDEX(mol->atoms, n[1]);
+ if (n[0] == n[1])
+ continue; /* Same atom */
+ if (exinfo != NULL) {
+ /* Look up exclusion table to exclude 1-2, 1-3, and 1-4 pairs */
+ for (i = exn1; i < exn2; i++) {
+ if (exlist[i] == n[1])
+ break;
+ }
+ if (i < exn2)
+ continue; /* Should be excluded */
+ }
+ if (MoleculeMeasureBond(mol, &r1, &(ap2->r)) < lim) {
+ /* Is this pair already registered? */
+ Int *ip;
+ for (i = 0, ip = pairs; i < npairs; i++, ip += 2) {
+ if ((ip[0] == n[0] && ip[1] == n[1]) || (ip[0] == n[1] && ip[1] == n[0]))
+ break;
+ }
+ if (i >= npairs) {
+ /* Not registered yet */
+ AssignArray(&pairs, &npairs, sizeof(Int) * 2, npairs, n);
+ }
+ }
+ }
+ }
+ 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;
}
/*
}
return INT2NUM(nbonds);
}
-
+
+#pragma mark ------ Cell and Symmetry ------
+
/*
* call-seq:
- * register_undo(script, *args)
+ * cell -> [a, b, c, alpha, beta, gamma [, sig_a, sig_b, sig_c, sig_alpha, sig_beta, sig_gamma]]
*
- * Register an undo operation with the current molecule.
+ * Returns the unit cell parameters. If cell is not set, returns nil.
*/
static VALUE
-s_Molecule_RegisterUndo(int argc, VALUE *argv, VALUE self)
+s_Molecule_Cell(VALUE self)
{
- Molecule *mol;
- VALUE script, args;
- MolAction *act;
+ Molecule *mol;
+ int i;
+ VALUE val;
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;
+ 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:
- * undo_enabled? -> bool
+ * 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)
*
- * Returns true if undo is enabled for this molecule; otherwise no.
+ * 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_UndoEnabled(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);
- if (MolActionCallback_isUndoRegistrationEnabled(mol))
- return Qtrue;
- else return Qfalse;
+ 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:
- * undo_enabled = bool
+ * box -> [avec, bvec, cvec, origin, flags]
*
- * Enable or disable undo.
+ * 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_SetUndoEnabled(VALUE self, VALUE val)
+s_Molecule_Box(VALUE self)
{
Molecule *mol;
+ VALUE v[5], val;
Data_Get_Struct(self, Molecule, mol);
- MolActionCallback_setUndoRegistrationEnabled(mol, (val != Qfalse && val != Qnil));
+ 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:
- * selection -> IntGroup
+ * 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
*
- * Returns the current selection.
+ * 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_Selection(VALUE self)
+s_Molecule_SetBox(VALUE self, VALUE aval)
{
Molecule *mol;
- IntGroup *ig;
- VALUE val;
+ 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);
- 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);
+ 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;
} else {
- val = IntGroup_Alloc(rb_cIntGroup);
+ 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;
}
- return val;
+ MolActionCreateAndPerform(mol, gMolActionSetBox, &(vv[0]), &(vv[1]), &(vv[2]), &origin, (flags[0] * 4 + flags[1] * 2 + flags[2]), convertCoordinates);
+ return self;
}
+/*
+ * call-seq:
+ * cell_periodicity -> [n1, n2, n3]
+ *
+ * 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_SetSelectionSub(VALUE self, VALUE val, int undoable)
+s_Molecule_CellPeriodicity(VALUE self)
{
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;
+ 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:
- * selection = IntGroup
+ * self.cell_periodicity = [n1, n2, n3] or Integer or nil
+ * set_cell_periodicity = [n1, n2, n3] or Integer or nil
*
- * 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.
+ * 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_SetSelection(VALUE self, VALUE val)
+s_Molecule_SetCellPeriodicity(VALUE self, VALUE arg)
{
- return s_Molecule_SetSelectionSub(self, val, 0);
+ 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;
}
/*
* call-seq:
- * set_undoable_selection(IntGroup)
+ * cell_flexibility -> bool
*
- * Set the current selection with undo registration. The right-hand operand may be nil.
- * This operation is undoable.
+ * Returns the unit cell is flexible or not
*/
static VALUE
-s_Molecule_SetUndoableSelection(VALUE self, VALUE val)
+s_Molecule_CellFlexibility(VALUE self)
{
- return s_Molecule_SetSelectionSub(self, val, 1);
+ 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:
- * hidden_atoms -> IntGroup
+ * self.cell_flexibility = bool
+ * set_cell_flexibility(bool)
*
- * Returns the currently hidden atoms.
+ * Change the unit cell is flexible or not
*/
static VALUE
-s_Molecule_HiddenAtoms(VALUE self)
+s_Molecule_SetCellFlexibility(VALUE self, VALUE arg)
{
- 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; */
+ 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_hidden_atoms(IntGroup)
- * self.hidden_atoms = IntGroup
+ * cell_transform -> Transform
*
- * Hide the specified atoms. This operation is _not_ undoable.
+ * Get the transform matrix that converts internal coordinates to cartesian coordinates.
+ * If cell is not defined, nil is returned.
*/
static VALUE
-s_Molecule_SetHiddenAtoms(VALUE self, VALUE val)
+s_Molecule_CellTransform(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;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol == NULL || mol->cell == NULL)
+ return Qnil;
+ return ValueFromTransform(&(mol->cell->tr));
+}
+
/*
- Molecule *mol;
+ * call-seq:
+ * symmetry -> Array of Transforms
+ * symmetries -> Array of Transforms
+ *
+ * Get the currently defined symmetry operations. If no symmetry operation is defined,
+ * returns an empty array.
+ */
+static VALUE
+s_Molecule_Symmetry(VALUE self)
+{
+ Molecule *mol;
+ VALUE val;
+ int i;
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);
+ 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; */
+ return val;
}
/*
* call-seq:
- * select_frame(index)
- * frame = index
+ * nsymmetries -> Integer
*
- * Select the specified frame. If successful, returns true, otherwise returns false.
+ * Get the number of currently defined symmetry operations.
*/
static VALUE
-s_Molecule_SelectFrame(VALUE self, VALUE val)
+s_Molecule_Nsymmetries(VALUE self)
{
Molecule *mol;
- int ival = NUM2INT(val);
Data_Get_Struct(self, Molecule, mol);
- ival = MoleculeSelectFrame(mol, ival, 1);
- if (ival >= 0)
- return Qtrue;
- else return Qfalse;
+ return INT2NUM(mol->nsyms);
}
/*
* call-seq:
- * frame -> Integer
+ * add_symmetry(Transform) -> Integer
*
- * Get the current frame.
+ * 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_Frame(VALUE self)
+s_Molecule_AddSymmetry(VALUE self, VALUE trans)
{
Molecule *mol;
+ Transform tr;
Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(mol->cframe);
+ TransformFromValue(trans, &tr);
+ MolActionCreateAndPerform(mol, gMolActionAddSymmetryOperation, &tr);
+ return INT2NUM(mol->nsyms);
}
/*
* call-seq:
- * nframes -> Integer
+ * remove_symmetry(count = nil) -> Integer
+ * remove_symmetries(count = nil) -> Integer
*
- * Get the number of frames.
+ * 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_Nframes(VALUE self)
+s_Molecule_RemoveSymmetry(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
+ VALUE cval;
+ int i, n;
Data_Get_Struct(self, Molecule, mol);
- return INT2NUM(MoleculeGetNumberOfFrames(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);
}
/*
* call-seq:
- * insert_frame(integer, coordinates = nil, cell_axes = nil) -> bool
- * insert_frames(intGroup = nil, coordinates = nil, cell_axes = nil) -> bool
+ * wrap_unit_cell(group) -> Vector3D
*
- * 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.
+ * 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_InsertFrames(int argc, VALUE *argv, VALUE self)
+s_Molecule_WrapUnitCell(VALUE self, VALUE gval)
{
- VALUE val, coords, cells;
Molecule *mol;
IntGroup *ig;
- 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, "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]);
- }
- }
- 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);
-}
-
-/*
- * call-seq:
- * remove_frames(IntGroup, wantCoordinates = false)
- *
- * 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_RemoveFrames(int argc, VALUE *argv, 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;
-}
-
-/*
- * call-seq:
- * each_frame {|n| ...}
- *
- * 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_EachFrame(VALUE self)
-{
- int i, cframe, nframes;
- Molecule *mol;
- 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));
- }
- MoleculeSelectFrame(mol, cframe, 1);
- }
- return self;
-}
-
-/*
- * call-seq:
- * get_coord_from_frame(index, group = nil)
- *
- * 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_GetCoordFromFrame(int argc, VALUE *argv, 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;
-}
-
-/*
- * call-seq:
- * reorder_frames(old_indices)
- *
- * 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_ReorderFrames(VALUE self, VALUE aval)
-{
- 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);
- return self;
-}
-
-/*
- * call-seq:
- * set_atom_attr(index, key, value)
- *
- * Set the atom attribute for the specified atom.
- * This operation is undoable.
- */
-static VALUE
-s_Molecule_SetAtomAttr(VALUE self, VALUE idx, VALUE key, VALUE val)
-{
- Molecule *mol;
- VALUE aref, oldval;
- 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;
-}
-
-/*
- * call-seq:
- * get_atom_attr(index, key)
- *
- * Get the atom attribute for the specified atom.
- */
-static VALUE
-s_Molecule_GetAtomAttr(VALUE self, VALUE idx, VALUE key)
-{
- return s_Molecule_SetAtomAttr(self, idx, key, Qundef);
-}
-
-/*
- * call-seq:
- * fragment(n1, *exatoms) -> IntGroup
- * fragment(group, *exatoms) -> IntGroup
- *
- * 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_Fragment(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- IntGroup *baseg, *ig, *exatoms;
- int n;
- volatile VALUE nval, exval;
+ Vector v, cv, dv;
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);
+ if (mol->cell == NULL)
+ rb_raise(rb_eMolbyError, "no unit cell is defined");
+ ig = s_Molecule_AtomGroupFromValue(self, gval);
+ s_Molecule_DoCenterOfMass(mol, &cv, ig);
+ TransformVec(&v, mol->cell->rtr, &cv);
+ if (mol->cell->flags[0])
+ v.x -= floor(v.x);
+ if (mol->cell->flags[1])
+ v.y -= floor(v.y);
+ if (mol->cell->flags[2])
+ v.z -= floor(v.z);
+ TransformVec(&dv, mol->cell->tr, &v);
+ VecDec(dv, cv);
+ MolActionCreateAndPerform(mol, gMolActionTranslateAtoms, &dv, ig);
IntGroupRelease(ig);
- return nval;
-}
-
-/*
- * call-seq:
- * fragments(exclude = nil)
- *
- * 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_Fragments(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- IntGroup *ag, *fg, *eg;
- VALUE gval, exval, retval;
- 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;
+ return ValueFromVector(&dv);
}
/*
* call-seq:
- * each_fragment(exclude = nil) {|group| ...}
+ * expand_by_symmetry(group, sym, dx=0, dy=0, dz=0, allow_overlap = false) -> Array
*
- * 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.
+ * 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_EachFragment(int argc, VALUE *argv, VALUE self)
+s_Molecule_ExpandBySymmetry(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- IntGroup *ag, *fg, *eg;
- VALUE gval, exval;
- 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);
- }
- IntGroupRelease(ag);
- if (eg != NULL)
- IntGroupRelease(eg);
- return self;
-}
-
-/*
- * call-seq:
- * detachable?(group) -> [n1, n2]
- *
- * 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_Detachable_P(VALUE self, VALUE gval)
-{
- Molecule *mol;
+ VALUE gval, sval, xval, yval, zval, rval, oval;
IntGroup *ig;
- int n1, n2;
- VALUE retval;
+ Int n[4], allow_overlap;
+ Int natoms;
+ Int nidx, *idx;
+
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 (MoleculeIsFragmentDetachable(mol, ig, &n1, &n2)) {
- retval = rb_ary_new3(2, INT2NUM(n1), INT2NUM(n2));
- } else retval = Qnil;
- IntGroupRelease(ig);
- return retval;
-}
-
-/*
- * call-seq:
- * bonds_on_border(group = selection) -> Array of Array of two Integers
- *
- * Returns an array of bonds that connect an atom in the group and an atom out
- * of the group. The first atom in the bond always belongs to the group. If no
- * such bonds are present, an empty array is returned.
- */
-static VALUE
-s_Molecule_BondsOnBorder(int argc, VALUE *argv, VALUE self)
-{
- Molecule *mol;
- IntGroup *ig, *bg;
- VALUE gval, retval;
- Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &gval);
- if (gval == Qnil) {
- ig = MoleculeGetSelection(mol);
- if (ig != NULL)
- IntGroupRetain(ig);
- } else {
- ig = s_Molecule_AtomGroupFromValue(self, gval);
- }
- retval = rb_ary_new();
- if (ig == NULL)
- return retval;
- bg = MoleculeSearchBondsAcrossAtomGroup(mol, ig);
- if (bg != NULL) {
- IntGroupIterator iter;
- Int i;
- IntGroupIteratorInit(bg, &iter);
- while ((i = IntGroupIteratorNext(&iter)) >= 0) {
- /* The atoms at the border */
- Int n1, n2;
- n1 = mol->bonds[i * 2];
- n2 = mol->bonds[i * 2 + 1];
- if (IntGroupLookupPoint(ig, n1) < 0) {
- int w = n1;
- n1 = n2;
- n2 = w;
- if (IntGroupLookupPoint(ig, n1) < 0)
- continue; /* Actually this is an internal error */
- }
- rb_ary_push(retval, rb_ary_new3(2, INT2NUM(n1), INT2NUM(n2)));
- }
- IntGroupIteratorRelease(&iter);
+ 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;
- IntGroup *ig;
- Vector vmin, vmax;
- int n;
- Atom *ap;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "01", &group);
- ig = (NIL_P(group) ? NULL : s_Molecule_AtomGroupFromValue(self, group));
- if (ig != NULL && IntGroupGetCount(ig) == 0)
- rb_raise(rb_eMolbyError, "atom group is empty");
- vmin.x = vmin.y = vmin.z = 1e30;
- vmax.x = vmax.y = vmax.z = -1e30;
- for (n = 0, ap = mol->atoms; n < mol->natoms; n++, ap = ATOM_NEXT(ap)) {
- Vector r;
- if (ig != NULL && IntGroupLookup(ig, n, NULL) == 0)
- continue;
- r = ap->r;
- if (r.x < vmin.x)
- vmin.x = r.x;
- if (r.y < vmin.y)
- vmin.y = r.y;
- if (r.z < vmin.z)
- vmin.z = r.z;
- if (r.x > vmax.x)
- vmax.x = r.x;
- if (r.y > vmax.y)
- vmax.y = r.y;
- if (r.z > vmax.z)
- vmax.z = r.z;
- }
- return rb_ary_new3(2, ValueFromVector(&vmin), ValueFromVector(&vmax));
+ return INT2NUM(mol->cframe);
}
-/* Get atom position or a vector */
-static void
-s_Molecule_GetVectorFromArg(Molecule *mol, VALUE val, Vector *vp)
+/*
+ * call-seq:
+ * nframes -> Integer
+ *
+ * Get the number of frames.
+ */
+static VALUE
+s_Molecule_Nframes(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;
+ 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;
+ 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, "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 {
- VectorFromValue(val, vp);
+ 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]);
+ }
}
+ 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:
- * measure_bond(n1, n2) -> Float
+ * create_frame(coordinates = nil) -> Integer
+ * create_frames(coordinates = nil) -> Integer
*
- * 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.
+ * Same as molecule.insert_frames(nil, coordinates).
*/
static VALUE
-s_Molecule_MeasureBond(VALUE self, VALUE nval1, VALUE nval2)
+s_Molecule_CreateFrames(int argc, VALUE *argv, VALUE self)
{
- Molecule *mol;
- Vector v1, v2;
- Data_Get_Struct(self, Molecule, mol);
- s_Molecule_GetVectorFromArg(mol, nval1, &v1);
- s_Molecule_GetVectorFromArg(mol, nval2, &v2);
- return rb_float_new(MoleculeMeasureBond(mol, &v1, &v2));
+ 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_angle(n1, n2, n3) -> Float
+ * remove_frames(IntGroup, wantCoordinates = false)
*
- * 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.
+ * 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_MeasureAngle(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3)
+s_Molecule_RemoveFrames(int argc, VALUE *argv, VALUE self)
{
+ VALUE val, flag;
+ VALUE retval;
Molecule *mol;
- Vector v1, v2, v3;
- Double d;
+ IntGroup *ig;
+ int count;
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);
+ 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_dihedral(n1, n2, n3, n4) -> Float
+ * each_frame {|n| ...}
*
- * 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.
+ * 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_MeasureDihedral(VALUE self, VALUE nval1, VALUE nval2, VALUE nval3, VALUE nval4)
+s_Molecule_EachFrame(VALUE self)
{
+ int i, cframe, nframes;
Molecule *mol;
- Vector v1, v2, v3, v4;
- Double d;
Data_Get_Struct(self, Molecule, mol);
- s_Molecule_GetVectorFromArg(mol, nval1, &v1);
- s_Molecule_GetVectorFromArg(mol, nval2, &v2);
- s_Molecule_GetVectorFromArg(mol, nval3, &v3);
- s_Molecule_GetVectorFromArg(mol, nval4, &v4);
- d = MoleculeMeasureDihedral(mol, &v1, &v2, &v3, &v4);
- if (isnan(d))
- return Qnil; /* Cannot define */
- else return rb_float_new(d);
+ 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:
- * expand_by_symmetry(group, sym, dx=0, dy=0, dz=0, allow_overlap = false) -> Array
+ * get_coord_from_frame(index, group = nil)
*
- * 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.
+ * 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_ExpandBySymmetry(int argc, VALUE *argv, VALUE self)
+s_Molecule_GetCoordFromFrame(int argc, VALUE *argv, VALUE self)
{
- Molecule *mol;
- VALUE gval, sval, xval, yval, zval, rval, oval;
+ Molecule *mol;
+ VALUE ival, gval, cval;
+ Int index, i, j, n, nn;
IntGroup *ig;
- Int n[4], allow_overlap;
- Int natoms;
- Int nidx, *idx;
-
+ IntGroupIterator iter;
+ Atom *ap;
+ Vector *vp;
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]));
+ 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 (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;
+ 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:
- * amend_by_symmetry(group = nil) -> IntGroup
+ * reorder_frames(old_indices)
*
- * Expand the specified part of the molecule by the given symmetry operation.
- * Returns an IntGroup containing the added atoms.
+ * 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_AmendBySymmetry(int argc, VALUE *argv, VALUE self)
+s_Molecule_ReorderFrames(VALUE self, VALUE aval)
{
- Molecule *mol;
- IntGroup *ig, *ig2;
- VALUE rval, gval;
+ Molecule *mol;
+ Int *ip, *ip2, i, n, nframes;
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;
+ 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);
+ return self;
}
+#pragma mark ------ Fragments ------
+
/*
* call-seq:
- * transform_for_symop(symop, is_cartesian = nil) -> Transform
+ * fragment(n1, *exatoms) -> IntGroup
+ * fragment(group, *exatoms) -> IntGroup
*
- * 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.
+ * 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_TransformForSymop(int argc, VALUE *argv, VALUE self)
+s_Molecule_Fragment(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE sval, fval;
- Symop symop;
- Transform tr;
+ IntGroup *baseg, *ig, *exatoms;
+ int n;
+ volatile VALUE nval, 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", &sval, &fval);
- if (rb_obj_is_kind_of(sval, rb_cNumeric)) {
- symop.sym = NUM2INT(rb_Integer(sval));
- symop.dx = symop.dy = symop.dz = 0;
+ 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 {
- 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]));
+ baseg = s_Molecule_AtomGroupFromValue(self, nval);
}
- 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);
+ 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:
- * symop_for_transform(transform, is_cartesian = nil) -> [sym, dx, dy, dz]
+ * fragments(exclude = nil)
*
- * 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.
+ * 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_SymopForTransform(int argc, VALUE *argv, VALUE self)
+s_Molecule_Fragments(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE tval, fval;
- Symop symop;
- Transform tr;
- int n;
+ 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", &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)
+ 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;
}
/*
* call-seq:
- * wrap_unit_cell(group) -> Vector3D
+ * each_fragment(exclude = nil) {|group| ...}
*
- * 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.
+ * 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_WrapUnitCell(VALUE self, VALUE gval)
+s_Molecule_EachFragment(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- IntGroup *ig;
- Vector v, cv, dv;
+ 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");
- ig = s_Molecule_AtomGroupFromValue(self, gval);
- s_Molecule_DoCenterOfMass(mol, &cv, ig);
- TransformVec(&v, mol->cell->rtr, &cv);
- if (mol->cell->flags[0])
- v.x -= floor(v.x);
- if (mol->cell->flags[1])
- v.y -= floor(v.y);
- if (mol->cell->flags[2])
- v.z -= floor(v.z);
- TransformVec(&dv, mol->cell->tr, &v);
- VecDec(dv, cv);
- MolActionCreateAndPerform(mol, gMolActionTranslateAtoms, &dv, ig);
- IntGroupRelease(ig);
- return ValueFromVector(&dv);
+ 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:
- * find_conflicts(limit[, group1[, group2 [, ignore_exclusion]]]) -> [[n1, n2], [n3, n4], ...]
+ * detachable?(group) -> [n1, n2]
*
- * 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.
+ * 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_FindConflicts(int argc, VALUE *argv, VALUE self)
+s_Molecule_Detachable_P(VALUE self, VALUE gval)
{
- 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;
+ int n1, n2;
+ VALUE 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;
+ ig = s_Molecule_AtomGroupFromValue(self, gval);
+ if (MoleculeIsFragmentDetachable(mol, ig, &n1, &n2)) {
+ retval = rb_ary_new3(2, INT2NUM(n1), INT2NUM(n2));
+ } else retval = Qnil;
+ IntGroupRelease(ig);
+ return retval;
+}
+
+/*
+ * call-seq:
+ * bonds_on_border(group = selection) -> Array of Array of two Integers
+ *
+ * Returns an array of bonds that connect an atom in the group and an atom out
+ * of the group. The first atom in the bond always belongs to the group. If no
+ * such bonds are present, an empty array is returned.
+ */
+static VALUE
+s_Molecule_BondsOnBorder(int argc, VALUE *argv, VALUE self)
+{
+ Molecule *mol;
+ IntGroup *ig, *bg;
+ VALUE gval, retval;
+ Data_Get_Struct(self, Molecule, mol);
+ rb_scan_args(argc, argv, "01", &gval);
+ if (gval == Qnil) {
+ ig = MoleculeGetSelection(mol);
+ if (ig != NULL)
+ IntGroupRetain(ig);
} else {
- 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
VALUE gval, rval, wval;
IntGroup *ig;
IntGroupIterator iter;
- int nn, errno, i, j, in, status;
+ int nn, errnum, i, j, in, status;
Vector *ref;
Double *weights, dval[3];
Transform tr;
if (gval == Qnil)
ig = IntGroupNewWithPoints(0, mol->natoms, -1);
else
- ig = IntGroupFromValue(gval);
+ ig = s_Molecule_AtomGroupFromValue(self, gval);
if (ig == NULL || (nn = IntGroupGetCount(ig)) == 0) {
IntGroupRelease(ig);
rb_raise(rb_eMolbyError, "atom group is not given correctly");
if (rb_obj_is_kind_of(rval, rb_cNumeric)) {
int fn = NUM2INT(rb_Integer(rval));
if (fn < 0 || fn >= MoleculeGetNumberOfFrames(mol)) {
- errno = 1;
+ errnum = 1;
status = fn;
goto err;
}
} else if (rb_obj_is_kind_of(rval, rb_cLAMatrix)) {
LAMatrix *m = LAMatrixFromValue(rval, NULL, 0, 0);
if (m->row * m->column < nn * 3) {
- errno = 2;
+ errnum = 2;
goto err;
}
for (i = 0; i < nn; i++) {
VALUE aval;
rval = rb_protect(rb_ary_to_ary, rval, &status);
if (status != 0) {
- errno = 3;
+ errnum = 3;
goto err;
}
if (RARRAY_LEN(rval) < nn) {
- errno = 2;
+ errnum = 2;
goto err;
}
if (rb_obj_is_kind_of((RARRAY_PTR(rval))[0], rb_cNumeric)) {
/* Array of 3*nn numbers */
if (RARRAY_LEN(rval) < nn * 3) {
- errno = 2;
+ errnum = 2;
goto err;
}
for (i = 0; i < nn; i++) {
for (j = 0; j < 3; j++) {
aval = rb_protect(rb_Float, (RARRAY_PTR(rval))[i * 3 + j], &status);
if (status != 0) {
- errno = 3;
+ errnum = 3;
goto err;
}
dval[j] = NUM2DBL(aval);
} else {
aval = rb_protect(rb_ary_to_ary, aval, &status);
if (status != 0) {
- errno = 3;
+ errnum = 3;
goto err;
}
if (RARRAY_LEN(aval) < 3) {
- errno = 4;
+ errnum = 4;
status = i;
goto err;
}
for (j = 0; j < 3; j++) {
VALUE aaval = rb_protect(rb_Float, (RARRAY_PTR(aval))[j], &status);
if (status != 0) {
- errno = 3;
+ errnum = 3;
goto err;
}
dval[j] = NUM2DBL(aaval);
} else {
wval = rb_protect(rb_ary_to_ary, wval, &status);
if (status != 0) {
- errno = 3;
+ errnum = 3;
goto err;
}
if (RARRAY_LEN(wval) < nn) {
- errno = 5;
+ errnum = 5;
goto err;
}
for (i = 0; i < nn; i++) {
VALUE wwval = rb_protect(rb_Float, (RARRAY_PTR(wval))[i], &status);
if (status != 0) {
- errno = 3;
+ errnum = 3;
goto err;
}
weights[i] = NUM2DBL(wwval);
}
dval[0] = s_Molecule_FitCoordinates_Sub(mol, ig, ref, weights, tr);
if (dval[0] < 0) {
- errno = 6;
+ errnum = 6;
goto err;
}
- errno = 0;
+ errnum = 0;
err:
IntGroupIteratorRelease(&iter);
free(ref);
free(weights);
- if (errno == 0) {
+ if (errnum == 0) {
return rb_ary_new3(2, ValueFromTransform(&tr), rb_float_new(dval[0]));
- } else if (errno == 1) {
+ } else if (errnum == 1) {
rb_raise(rb_eMolbyError, "frame index (%d) is out of range", status);
- } else if (errno == 2) {
+ } else if (errnum == 2) {
rb_raise(rb_eMolbyError, "insufficient number of reference coordinates");
- } else if (errno == 3) {
+ } else if (errnum == 3) {
rb_jump_tag(status);
- } else if (errno == 4) {
+ } else if (errnum == 4) {
rb_raise(rb_eMolbyError, "less than 3 elements for index %d of reference coordinates", status);
- } else if (errno == 5) {
+ } else if (errnum == 5) {
rb_raise(rb_eMolbyError, "insufficient number of weight values");
- } else if (errno == 6) {
+ } else if (errnum == 6) {
rb_raise(rb_eMolbyError, "matrix calculation failed during coordinate fitting");
}
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:
- * create_graphic(kind, color, points, fill = nil) -> integer
+ * export_graphic(fname, scale = 1.0, bg_color = -1, width = 0, height = 0)
+ *
+ * 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.
+ * If either width or height is not specified, then the screen width/height is used instead.
+ */
+static VALUE
+s_Molecule_ExportGraphic(int argc, VALUE *argv, VALUE self)
+{
+ Molecule *mol;
+ VALUE fval, sval, bval, wval, hval;
+ char *fname;
+ float scale;
+ int bg_color, width, height;
+ 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, "14", &fval, &sval, &bval, &wval, &hval);
+ 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 (wval == Qnil)
+ width = 0;
+ else width = NUM2INT(rb_Integer(wval));
+ if (hval == Qnil)
+ height = 0;
+ else height = NUM2INT(rb_Integer(hval));
+ if (MainViewCallback_exportGraphic(mol->mview, fname, scale, bg_color, width, height) == 0)
+ return fval;
+ else return Qnil;
+}
+
+#pragma mark ------ Graphics ------
+
+static void
+s_CalculateGraphicNormals(MainViewGraphic *gp)
+{
+ int i;
+ Vector v1, v2, v3;
+ if (gp == NULL || gp->npoints < 3)
+ return;
+ AssignArray(&gp->normals, &gp->nnormals, sizeof(GLfloat) * 3, gp->npoints - 1, NULL);
+ v1.x = gp->points[3] - gp->points[0];
+ v1.y = gp->points[4] - gp->points[1];
+ v1.z = gp->points[5] - gp->points[2];
+ /* nv[i] = (v[i-1]-v[0]).cross(v[i]-v[0]) (i=2..n-1) */
+ for (i = 2; i < gp->npoints; i++) {
+ v2.x = gp->points[i * 3] - gp->points[0];
+ v2.y = gp->points[i * 3 + 1] - gp->points[1];
+ v2.z = gp->points[i * 3 + 2] - gp->points[2];
+ VecCross(v3, v1, v2);
+ NormalizeVec(&v3, &v3);
+ gp->normals[i * 3] = v3.x;
+ gp->normals[i * 3 + 1] = v3.y;
+ gp->normals[i * 3 + 2] = v3.z;
+ v1 = v2;
+ }
+ /* normals[0] = average of all nv[i] (i=2..n-1) */
+ VecZero(v1);
+ for (i = 2; i < gp->npoints; i++) {
+ v1.x += gp->normals[i * 3];
+ v1.y += gp->normals[i * 3 + 1];
+ v1.z += gp->normals[i * 3 + 2];
+ }
+ NormalizeVec(&v1, &v1);
+ gp->normals[0] = v1.x;
+ gp->normals[1] = v1.y;
+ gp->normals[2] = v1.z;
+ /* normals[1] = nv[2].normalize */
+ v2.x = gp->normals[6];
+ v2.y = gp->normals[7];
+ v2.z = gp->normals[8];
+ NormalizeVec(&v1, &v2);
+ gp->normals[3] = v1.x;
+ gp->normals[4] = v1.y;
+ gp->normals[5] = v1.z;
+ /* normals[i] = (nv[i] + nv[i+1]).normalize (i=2..n-2) */
+ for (i = 2; i < gp->npoints; i++) {
+ if (i == gp->npoints - 1)
+ VecZero(v3);
+ else {
+ v3.x = gp->normals[i * 3 + 3];
+ v3.y = gp->normals[i * 3 + 4];
+ v3.z = gp->normals[i * 3 + 5];
+ }
+ VecInc(v2, v3);
+ NormalizeVec(&v1, &v2);
+ gp->normals[i * 3] = v1.x;
+ gp->normals[i * 3 + 1] = v1.y;
+ gp->normals[i * 3 + 2] = v1.z;
+ v2 = v3;
+ }
+}
+
+/*
+ * call-seq:
+ * insert_graphic(index, kind, color, points, fill = nil) -> integer
*
- * Create a new graphic object.
+ * Create a new graphic object and insert at the given graphic index (if -1, then append at the last).
* kind: a symbol representing the kind of the graphic. :line, :poly, :cylinder, :cone, :ellipsoid
* color: an array of 3 (rgb) or 4 (rgba) floating numbers
* points: an array of Vectors
*
*/
static VALUE
-s_Molecule_CreateGraphic(int argc, VALUE *argv, VALUE self)
+s_Molecule_InsertGraphic(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
MainViewGraphic g;
- int i, n, ni;
+ int i, n, ni, idx;
const char *p;
- VALUE kval, cval, pval, fval;
+ VALUE kval, cval, pval, fval, ival;
Data_Get_Struct(self, Molecule, mol);
if (mol->mview == NULL)
rb_raise(rb_eMolbyError, "this molecule has no associated graphic view");
- rb_scan_args(argc, argv, "31", &kval, &cval, &pval, &fval);
- kval = rb_obj_as_string(kval);
+ rb_scan_args(argc, argv, "41", &ival, &kval, &cval, &pval, &fval);
+ idx = NUM2INT(rb_Integer(ival));
+ if (idx == -1)
+ idx = mol->mview->ngraphics;
+ else if (idx < 0 || idx > mol->mview->ngraphics)
+ rb_raise(rb_eMolbyError, "the graphic index (%d) out of range", idx);
memset(&g, 0, sizeof(g));
g.visible = 1;
- p = RSTRING_PTR(kval);
- if (strcmp(p, "line") == 0)
- g.kind = kMainViewGraphicLine;
- else if (strcmp(p, "poly") == 0)
- g.kind = kMainViewGraphicPoly;
- else if (strcmp(p, "cylinder") == 0)
- g.kind = kMainViewGraphicCylinder;
- else if (strcmp(p, "cone") == 0)
- g.kind = kMainViewGraphicCone;
- else if (strcmp(p, "ellipsoid") == 0)
- g.kind = kMainViewGraphicEllipsoid;
- else rb_raise(rb_eMolbyError, "unknown graphic object type: %s", p);
+ if (rb_obj_is_kind_of(kval, rb_cInteger)) {
+ g.kind = NUM2INT(kval); /* Direct assign (for undo registration) */
+ } else {
+ kval = rb_obj_as_string(kval);
+ p = StringValuePtr(kval);
+ if (strcmp(p, "line") == 0)
+ g.kind = kMainViewGraphicLine;
+ else if (strcmp(p, "poly") == 0)
+ g.kind = kMainViewGraphicPoly;
+ else if (strcmp(p, "cylinder") == 0)
+ g.kind = kMainViewGraphicCylinder;
+ else if (strcmp(p, "cone") == 0)
+ g.kind = kMainViewGraphicCone;
+ else if (strcmp(p, "ellipsoid") == 0)
+ g.kind = kMainViewGraphicEllipsoid;
+ else rb_raise(rb_eMolbyError, "unknown graphic object type: %s", p);
+ }
g.closed = (RTEST(fval) ? 1 : 0);
cval = rb_ary_to_ary(cval);
n = RARRAY_LEN(cval);
NewArray(&g.points, &g.npoints, sizeof(GLfloat) * 3, n);
for (i = 0; i < n; i++) {
Vector v;
+ VALUE rval = RARRAY_PTR(pval)[i];
if (i == ni) {
- v.x = NUM2DBL(rb_Float(RARRAY_PTR(pval)[i]));
+ if (rb_obj_is_kind_of(rval, rb_cVector3D)) {
+ /* The float argument can also be given as a vector (for simplify undo registration) */
+ VectorFromValue(rval, &v);
+ } else {
+ v.x = NUM2DBL(rb_Float(rval));
+ }
v.y = v.z = 0;
} else {
- VectorFromValue(RARRAY_PTR(pval)[i], &v);
+ VectorFromValue(rval, &v);
}
g.points[i * 3] = v.x;
g.points[i * 3 + 1] = v.y;
g.points[6] = g.points[8] = g.points[9] = g.points[10] = 0;
g.points[7] = g.points[11] = g.points[3];
}
- MainView_insertGraphic(mol->mview, -1, &g);
- return INT2NUM(mol->mview->ngraphics - 1);
+ if (g.kind == kMainViewGraphicPoly) {
+ /* Calculate normals */
+ s_CalculateGraphicNormals(&g);
+ }
+ MainView_insertGraphic(mol->mview, idx, &g);
+
+ {
+ /* Register undo */
+ MolAction *act;
+ act = MolActionNew(SCRIPT_ACTION("i"), "remove_graphic", idx);
+ MolActionCallback_registerUndo(mol, act);
+ MolActionRelease(act);
+ }
+
+ return INT2NUM(idx);
+}
+
+/*
+ * call-seq:
+ * create_graphic(kind, color, points, fill = nil) -> integer
+ *
+ * Create a new graphic object. The arguments are similar as insert_graphic.
+ */
+static VALUE
+s_Molecule_CreateGraphic(int argc, VALUE *argv, VALUE self)
+{
+ VALUE args[5];
+ rb_scan_args(argc, argv, "31", args + 1, args + 2, args + 3, args + 4);
+ args[0] = INT2NUM(-1);
+ return s_Molecule_InsertGraphic(argc + 1, args, self);
}
/*
i = NUM2INT(rb_Integer(ival));
if (i < 0 || i >= mol->mview->ngraphics)
rb_raise(rb_eArgError, "graphic index is out of range");
+ {
+ /* Prepare data for undo */
+ MainViewGraphic *gp;
+ Vector *vp;
+ MolAction *act;
+ double col[4];
+ int n;
+ gp = mol->mview->graphics + i;
+ vp = (Vector *)malloc(sizeof(Vector) * gp->npoints);
+ for (n = 0; n < gp->npoints; n++) {
+ vp[n].x = gp->points[n * 3];
+ vp[n].y = gp->points[n * 3 + 1];
+ vp[n].z = gp->points[n * 3 + 2];
+ }
+ col[0] = gp->rgba[0];
+ col[1] = gp->rgba[1];
+ col[2] = gp->rgba[2];
+ col[3] = gp->rgba[3];
+ if (gp->visible == 0) {
+ act = MolActionNew(SCRIPT_ACTION("i"), "hide_graphic", i);
+ MolActionCallback_registerUndo(mol, act);
+ MolActionRelease(act);
+ }
+ act = MolActionNew(SCRIPT_ACTION("iiDVb"), "insert_graphic", i, gp->kind, 4, col, gp->npoints, vp, gp->closed);
+ MolActionCallback_registerUndo(mol, act);
+ free(vp);
+ MolActionRelease(act);
+ }
MainView_removeGraphic(mol->mview, i);
return ival;
}
rb_raise(rb_eMolbyError, "this molecule has no associated graphic view");
return INT2NUM(mol->mview->ngraphics);
}
-
+
+/*
+ * call-seq:
+ * get_graphic_point(graphic_index, point_index) -> value
+ * get_graphic_points(graphic_index) -> values
+ *
+ * Get the point_index-th control point of graphic_index-th graphic object.
+ * Get an array of all control points with the given values.
+ *
+ */
+static VALUE
+s_Molecule_GetGraphicPoint(int argc, VALUE *argv, VALUE self)
+{
+ MainViewGraphic *gp;
+ Molecule *mol;
+ int index, pindex;
+ Vector v;
+ VALUE gval, pval;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->mview == NULL)
+ rb_raise(rb_eMolbyError, "this molecule has no associated graphic view");
+ rb_scan_args(argc, argv, "11", &gval, &pval);
+ index = NUM2INT(rb_Integer(gval));
+ if (index < 0 || index >= mol->mview->ngraphics)
+ rb_raise(rb_eArgError, "the graphic index is out of range");
+ gp = mol->mview->graphics + index;
+ if (pval != Qnil) {
+ pindex = NUM2INT(rb_Integer(pval));
+ if (pindex < 0 || pindex >= gp->npoints)
+ rb_raise(rb_eArgError, "the point index is out of range");
+ v.x = gp->points[pindex * 3];
+ v.y = gp->points[pindex * 3 + 1];
+ v.z = gp->points[pindex * 3 + 2];
+ if ((gp->kind == kMainViewGraphicCylinder || gp->kind == kMainViewGraphicCone) && pindex == 2) {
+ return rb_float_new(v.x);
+ } else {
+ return ValueFromVector(&v);
+ }
+ } else {
+ pval = rb_ary_new();
+ for (pindex = 0; pindex < gp->npoints; pindex++) {
+ v.x = gp->points[pindex * 3];
+ v.y = gp->points[pindex * 3 + 1];
+ v.z = gp->points[pindex * 3 + 2];
+ rb_ary_push(pval, ValueFromVector(&v));
+ }
+ return pval;
+ }
+}
+
+/*
+ * call-seq:
+ * set_graphic_point(graphic_index, point_index, new_value) -> new_value
+ * set_graphic_points(graphic_index, new_values) -> new_values
+ *
+ * Change the point_index-th control point of graphic_index-th graphic object.
+ * Replace the control points with the given values.
+ *
+ */
+static VALUE
+s_Molecule_SetGraphicPoint(int argc, VALUE *argv, VALUE self)
+{
+ MainViewGraphic *gp;
+ Molecule *mol;
+ int index, pindex;
+ Vector v, v0;
+ VALUE gval, pval, nval;
+ MolAction *act;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->mview == NULL)
+ rb_raise(rb_eMolbyError, "this molecule has no associated graphic view");
+ rb_scan_args(argc, argv, "21", &gval, &pval, &nval);
+ index = NUM2INT(rb_Integer(gval));
+ if (index < 0 || index >= mol->mview->ngraphics)
+ rb_raise(rb_eArgError, "the graphic index is out of range");
+ gp = mol->mview->graphics + index;
+ if (nval != Qnil) {
+ pindex = NUM2INT(rb_Integer(pval));
+ if (pindex < 0 || pindex >= gp->npoints)
+ rb_raise(rb_eArgError, "the point index is out of range");
+ v0.x = gp->points[pindex * 3];
+ v0.y = gp->points[pindex * 3 + 1];
+ v0.z = gp->points[pindex * 3 + 2];
+ if (rb_obj_is_kind_of(nval, rb_cNumeric)) {
+ if ((gp->kind == kMainViewGraphicCylinder || gp->kind == kMainViewGraphicCone) && pindex == 2) {
+ v.x = NUM2DBL(rb_Float(nval));
+ v.y = v.z = 0;
+ } else if (gp->kind == kMainViewGraphicEllipsoid && pindex == 1) {
+ v.x = NUM2DBL(rb_Float(nval));
+ v.y = v.z = 0;
+ gp->points[7] = gp->points[11] = v.x;
+ gp->points[6] = gp->points[8] = gp->points[9] = gp->points[10] = 0;
+ } else rb_raise(rb_eArgError, "the argument must be an array-like object");
+ } else {
+ if (nval == Qnil) {
+ v.x = kInvalidFloat;
+ v.y = v.z = 0.0;
+ } else VectorFromValue(nval, &v);
+ }
+ gp->points[pindex * 3] = v.x;
+ gp->points[pindex * 3 + 1] = v.y;
+ gp->points[pindex * 3 + 2] = v.z;
+ act = MolActionNew(SCRIPT_ACTION("iiv"), "set_graphic_point", index, pindex, &v0);
+ } else {
+ VALUE aval;
+ int len;
+ Vector *vp = (Vector *)malloc(sizeof(Vector) * gp->npoints);
+ for (pindex = 0; pindex < gp->npoints; pindex++) {
+ vp[pindex].x = gp->points[pindex * 3];
+ vp[pindex].y = gp->points[pindex * 3 + 1];
+ vp[pindex].z = gp->points[pindex * 3 + 2];
+ }
+ act = MolActionNew(SCRIPT_ACTION("iV"), "set_graphic_points", index, gp->npoints, vp);
+ free(vp);
+ pval = rb_ary_to_ary(pval);
+ len = RARRAY_LEN(pval);
+ if (gp->npoints < len) {
+ gp->points = (GLfloat *)realloc(gp->points, sizeof(GLfloat) * 3 * len);
+ gp->npoints = len;
+ } else if (gp->npoints > len) {
+ int len2 = 3;
+ switch (gp->kind) {
+ case kMainViewGraphicLine: len2 = 2; break;
+ case kMainViewGraphicPoly: len2 = 3; break;
+ case kMainViewGraphicCylinder: len2 = 3; break;
+ case kMainViewGraphicCone: len2 = 3; break;
+ case kMainViewGraphicEllipsoid: len2 = 4; break;
+ }
+ if (len2 < len)
+ len2 = len;
+ gp->npoints = len2;
+ }
+ for (pindex = 0; pindex < len && pindex < gp->npoints; pindex++) {
+ aval = RARRAY_PTR(pval)[pindex];
+ if ((gp->kind == kMainViewGraphicCylinder || gp->kind == kMainViewGraphicCone) && pindex == 2) {
+ v.x = NUM2DBL(rb_Float(aval));
+ v.y = v.z = 0;
+ } else if (gp->kind == kMainViewGraphicEllipsoid && pindex == 1 && len == 2) {
+ v.x = NUM2DBL(rb_Float(aval));
+ v.y = v.z = 0;
+ gp->points[7] = gp->points[11] = v.x;
+ gp->points[6] = gp->points[8] = gp->points[9] = gp->points[10] = 0;
+ break;
+ } else VectorFromValue(aval, &v);
+ gp->points[pindex * 3] = v.x;
+ gp->points[pindex * 3 + 1] = v.y;
+ gp->points[pindex * 3 + 2] = v.z;
+ }
+ }
+ if (gp->kind == kMainViewGraphicPoly) {
+ /* Calculate normals */
+ s_CalculateGraphicNormals(gp);
+ }
+ MolActionCallback_registerUndo(mol, act);
+ MolActionRelease(act);
+ MoleculeCallback_notifyModification(mol, 0);
+ return nval;
+}
+
/*
* call-seq:
- * set_graphic_point(graphic_index, point_index, new_value) -> new_value
+ * get_graphic_color(graphic_index) -> value
*
- * Change the point_index-th control point of graphic_index-th graphic object
- *
+ * Get the color of graphic_index-th graphic object
*/
static VALUE
-s_Molecule_SetGraphicPoint(VALUE self, VALUE gval, VALUE pval, VALUE nval)
+s_Molecule_GetGraphicColor(VALUE self, VALUE gval)
{
MainViewGraphic *gp;
Molecule *mol;
int index;
- Vector v;
Data_Get_Struct(self, Molecule, mol);
if (mol->mview == NULL)
rb_raise(rb_eMolbyError, "this molecule has no associated graphic view");
if (index < 0 || index >= mol->mview->ngraphics)
rb_raise(rb_eArgError, "the graphic index is out of range");
gp = mol->mview->graphics + index;
- index = NUM2INT(rb_Integer(pval));
- if (index < 0 || index >= gp->npoints)
- rb_raise(rb_eArgError, "the point index is out of range");
- if (rb_obj_is_kind_of(nval, rb_cNumeric)) {
- if ((gp->kind == kMainViewGraphicCylinder || gp->kind == kMainViewGraphicCone) && index == 2) {
- v.x = NUM2DBL(rb_Float(nval));
- v.y = v.z = 0;
- } else if (gp->kind == kMainViewGraphicEllipsoid && index == 1) {
- gp->points[3] = gp->points[7] = gp->points[11] = NUM2DBL(rb_Float(nval));
- gp->points[4] = gp->points[5] = gp->points[6] = gp->points[8] = gp->points[9] = gp->points[10] = 0;
- return nval;
- } else rb_raise(rb_eArgError, "the argument must be an array-like object");
- } else {
- if (nval == Qnil) {
- v.x = kInvalidFloat;
- v.y = v.z = 0.0;
- } else VectorFromValue(nval, &v);
- }
- gp->points[index * 3] = v.x;
- gp->points[index * 3 + 1] = v.y;
- gp->points[index * 3 + 2] = v.z;
- MoleculeCallback_notifyModification(mol, 0);
- return nval;
+ return rb_ary_new3(4, rb_float_new(gp->rgba[0]), rb_float_new(gp->rgba[1]), rb_float_new(gp->rgba[2]), rb_float_new(gp->rgba[3]));
}
/*
{
MainViewGraphic *gp;
Molecule *mol;
- int index, n;
+ MolAction *act;
+ double c[4];
+ int index, i, n;
Data_Get_Struct(self, Molecule, mol);
if (mol->mview == NULL)
rb_raise(rb_eMolbyError, "this molecule has no associated graphic view");
if (index < 0 || index >= mol->mview->ngraphics)
rb_raise(rb_eArgError, "the graphic index is out of range");
gp = mol->mview->graphics + index;
+ for (i = 0; i < 4; i++)
+ c[i] = gp->rgba[i];
cval = rb_ary_to_ary(cval);
n = RARRAY_LEN(cval);
if (n != 3 && n != 4)
rb_raise(rb_eArgError, "the color argument must have 3 or 4 numbers");
- for (index = 0; index < n; index++) {
- gp->rgba[index] = NUM2DBL(rb_Float(RARRAY_PTR(cval)[index]));
+
+ for (i = 0; i < n; i++) {
+ gp->rgba[i] = NUM2DBL(rb_Float(RARRAY_PTR(cval)[i]));
}
if (n == 3)
gp->rgba[3] = 1.0;
+ act = MolActionNew(SCRIPT_ACTION("iD"), "set_graphic_color", index, 4, c);
+ MolActionCallback_registerUndo(mol, act);
+ MolActionRelease(act);
MoleculeCallback_notifyModification(mol, 0);
return cval;
}
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
/* Set up parameters */
mono = NUM2INT(rb_Integer(mval));
- if (mono <= 0 || mono > mol->bset->ncomps)
- rb_raise(rb_eMolbyError, "The MO number (%d) is out of range (should be 1..%d)", mono, mol->bset->ncomps);
+ if (mono < 0 || mono > mol->bset->ncomps)
+ rb_raise(rb_eMolbyError, "The MO number (%d) is out of range (should be 1..%d, or 0 as 'arbitrary vector')", mono, mol->bset->ncomps);
if (RTEST(bval)) {
if (mol->bset->rflag != 0)
rb_raise(rb_eMolbyError, "Beta MO is requested but not present");
/*
* call-seq:
+ * clear_surface
+ *
+ * Clear the MO surface if present.
+ */
+static VALUE
+s_Molecule_ClearSurface(VALUE self)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->mcube != NULL)
+ MoleculeClearMCube(mol, 0, 0, 0, NULL, 0.0, 0.0, 0.0);
+ return self;
+}
+
+/*
+ * call-seq:
+ * hide_surface
+ *
+ * Hide the MO surface if present.
+ */
+static VALUE
+s_Molecule_HideSurface(VALUE self)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->mcube != NULL) {
+ mol->mcube->hidden = 1;
+ MoleculeCallback_notifyModification(mol, 0);
+ }
+ return self;
+}
+
+/*
+ * call-seq:
+ * show_surface
+ *
+ * Show the MO surface if present.
+ */
+static VALUE
+s_Molecule_ShowSurface(VALUE self)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->mcube != NULL) {
+ mol->mcube->hidden = 0;
+ MoleculeCallback_notifyModification(mol, 0);
+ }
+ return self;
+}
+
+/*
+ * call-seq:
* create_surface(mo, attr = nil)
*
* Create a MO surface. The argument mo is the MO index (1-based); if mo is negative,
/*
* call-seq:
- * add_gaussian_orbital_shell(atom_index, sym, no_of_primitives)
+ * add_gaussian_orbital_shell(atom_index, sym, no_of_primitives[, additional_exponent])
*
* To be used internally. Add a gaussian orbital shell with the atom index, symmetry code,
- * and the number of primitives. Symmetry code: 0, S-type; 1, P-type; -1, SP-type; 2, D-type;
- * -2, D5-type.
+ * and the number of primitives.
+ * Additional exponent is for JANPA only; implements an additinal r^N component that
+ * appears in cartesian->spherical conversion.
+ * Symmetry code: 0, S-type; 1, P-type; -1, SP-type; 2, D-type; -2, D5-type;
+ * 3, F-type; -3, F7-type; 4, G-type; -4, G9-type.
+ * Or: "s", S-type; "p", P-type; "sp", SP-type; "d", D-type; "d5", D5-type;
+ * "f", F-type; "f7", F7-type; "g", G-type; "g9", G9-type
*/
static VALUE
-s_Molecule_AddGaussianOrbitalShell(VALUE self, VALUE aval, VALUE symval, VALUE npval)
+s_Molecule_AddGaussianOrbitalShell(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- int sym, nprims, a_idx, n;
- Data_Get_Struct(self, Molecule, mol);
+ int sym, nprims, a_idx, n, add_exp;
+ VALUE aval, symval, npval, addval;
+ Data_Get_Struct(self, Molecule, mol);
+ rb_scan_args(argc, argv, "31", &aval, &symval, &npval, &addval);
+ if (rb_obj_is_kind_of(symval, rb_cString)) {
+ const char *p = StringValuePtr(symval);
+ if (strcasecmp(p, "s") == 0)
+ sym = 0;
+ else if (strcasecmp(p, "p") == 0)
+ sym = 1;
+ else if (strcasecmp(p, "sp") == 0)
+ sym = -1;
+ else if (strcasecmp(p, "d") == 0)
+ sym = 2;
+ else if (strcasecmp(p, "d5") == 0)
+ sym = -2;
+ else if (strcasecmp(p, "f") == 0)
+ sym = 3;
+ else if (strcasecmp(p, "f7") == 0)
+ sym = -3;
+ else if (strcasecmp(p, "g") == 0)
+ sym = 4;
+ else if (strcasecmp(p, "g9") == 0)
+ sym = -4;
+ else
+ rb_raise(rb_eArgError, "Unknown orbital type '%s'", p);
+ } else {
+ sym = NUM2INT(rb_Integer(symval));
+ }
a_idx = NUM2INT(rb_Integer(aval));
- sym = NUM2INT(rb_Integer(symval));
nprims = NUM2INT(rb_Integer(npval));
- n = MoleculeAddGaussianOrbitalShell(mol, a_idx, sym, nprims);
+ if (addval != Qnil)
+ add_exp = NUM2INT(rb_Integer(addval));
+ else add_exp = 0;
+ n = MoleculeAddGaussianOrbitalShell(mol, a_idx, sym, nprims, add_exp);
if (n == -1)
rb_raise(rb_eMolbyError, "Molecule is emptry");
else if (n == -2)
n = MoleculeGetGaussianComponentInfo(mol, c, &atom_idx, label, &shell_idx);
if (n != 0)
rb_raise(rb_eMolbyError, "Cannot get the shell info for component index (%d)", c);
- return rb_ary_new3(3, INT2NUM(atom_idx), INT2NUM(shell_idx), rb_str_new2(label));
+ return rb_ary_new3(3, INT2NUM(atom_idx), INT2NUM(shell_idx), Ruby_NewEncodedStringValue2(label));
}
/*
return Qnil;
if (kval == sTypeSym) {
switch (mol->bset->rflag) {
- case 0: return rb_str_new2("UHF");
- case 1: return rb_str_new2("RHF");
- case 2: return rb_str_new2("ROHF");
+ case 0: return Ruby_NewEncodedStringValue2("UHF");
+ case 1: return Ruby_NewEncodedStringValue2("RHF");
+ case 2: return Ruby_NewEncodedStringValue2("ROHF");
default: return rb_str_to_str(INT2NUM(mol->bset->rflag));
}
} else if (kval == sAlphaSym) {
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;
-}
-#endif
+#pragma mark ------ Molecular Topology ------
/*
* call-seq:
return Qnil;
rb_scan_args(argc, argv, "01", &val);
if (val != Qnil)
- ig = IntGroupFromValue(val);
+ ig = s_Molecule_AtomGroupFromValue(self, val);
else ig = NULL;
result = MoleculeSearchEquivalentAtoms(mol, ig);
if (result == NULL)
gval = *argv++;
argc -= 2;
Data_Get_Struct(self, Molecule, mol);
- ig = IntGroupFromValue(gval);
+ if (gval == Qnil)
+ ig = NULL;
+ else
+ ig = s_Molecule_AtomGroupFromValue(self, gval);
+ if (ig == NULL || IntGroupGetCount(ig) == 0)
+ rb_raise(rb_eMolbyError, "atom group is not given correctly");
memset(&a, 0, sizeof(a));
memset(&an, 0, sizeof(an));
strncpy(a.aname, StringValuePtr(nval), 4);
return Data_Wrap_Struct(rb_cAtomRef, 0, (void (*)(void *))AtomRefRelease, aref);
}
+#pragma mark ------ Molecular Properties ------
+
/*
* call-seq:
* set_property(name, value[, index]) -> value
Data_Get_Struct(self, Molecule, mol);
rval = rb_ary_new();
for (i = mol->nmolprops - 1; i >= 0; i--) {
- nval = rb_str_new2(mol->molprops[i].propname);
+ nval = Ruby_NewEncodedStringValue2(mol->molprops[i].propname);
rb_ary_store(rval, i, nval);
}
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:
for (idx = 0; (mol = MoleculeCallback_moleculeAtIndex(idx)) != NULL; idx++) {
VALUE name;
MoleculeCallback_displayName(mol, buf, sizeof buf);
- name = rb_str_new2(buf);
+ name = Ruby_NewEncodedStringValue2(buf);
if (rb_reg_match(val, name) != Qnil && --k == 0)
break;
}
int i;
VALUE ary;
i = 0;
- ary = rb_ary_new();
- while ((mol = MoleculeCallback_moleculeAtIndex(i)) != NULL) {
- rb_ary_push(ary, ValueFromMolecule(mol));
- i++;
- }
- return ary;
-}
-
-/*
- * call-seq:
- * ordered_list -> array of Molecules
- *
- * Get the list of molecules associated to the documents, in the order of front-to-back
- * ordering of the associated window. If no document is open, returns an empry array.
- */
-static VALUE
-s_Molecule_OrderedList(VALUE klass)
-{
- Molecule *mol;
- int i;
- VALUE ary;
- i = 0;
- ary = rb_ary_new();
- while ((mol = MoleculeCallback_moleculeAtOrderedIndex(i)) != NULL) {
- rb_ary_push(ary, ValueFromMolecule(mol));
- i++;
- }
- 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));
+ ary = rb_ary_new();
+ while ((mol = MoleculeCallback_moleculeAtIndex(i)) != NULL) {
+ rb_ary_push(ary, ValueFromMolecule(mol));
+ i++;
}
- return sval;
+ return ary;
}
/*
* call-seq:
- * self == Molecule -> boolean
+ * ordered_list -> array of Molecules
*
- * True if the two arguments point to the same molecule.
+ * Get the list of molecules associated to the documents, in the order of front-to-back
+ * ordering of the associated window. If no document is open, returns an empry array.
*/
static VALUE
-s_Molecule_Equal(VALUE self, VALUE val)
+s_Molecule_OrderedList(VALUE klass)
{
- 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;
+ Molecule *mol;
+ int i;
+ VALUE ary;
+ i = 0;
+ ary = rb_ary_new();
+ while ((mol = MoleculeCallback_moleculeAtOrderedIndex(i)) != NULL) {
+ rb_ary_push(ary, ValueFromMolecule(mol));
+ i++;
+ }
+ return ary;
}
+#pragma mark ------ Call Subprocess ------
+
/* The callback functions for call_subprocess_async */
static int
s_Molecule_CallSubProcessAsync_EndCallback(Molecule *mol, int status)
* call_subprocess_async(cmd [, end_callback [, timer_callback [, standard_output_file [, error_output_file]]]])
*
* Call subprocess asynchronically.
+ * cmd is either a single string of an array of string. If it is a single string, then
+ * it will be given to wxExecute as a single argument. In this case, the string can be
+ * split into arguments by whitespace. If this behavior is not intended, then use an array
+ * containing a single string.
* If end_callback is given, it will be called (with two arguments self and termination status)
* when the subprocess terminated.
* If timer_callback is given, it will be called (also with two arguments, self and timer count).
VALUE cmd, end_proc, timer_proc, stdout_val, stderr_val;
Molecule *mol;
char *sout, *serr;
+ const char **cmdargv;
int n;
FILE *fpout, *fperr;
rb_scan_args(argc, argv, "14", &cmd, &end_proc, &timer_proc, &stdout_val, &stderr_val);
/* Register procs as instance variables */
rb_ivar_set(self, rb_intern("end_proc"), end_proc);
rb_ivar_set(self, rb_intern("timer_proc"), timer_proc);
- n = MoleculeCallback_callSubProcessAsync(mol, StringValuePtr(cmd), s_Molecule_CallSubProcessAsync_EndCallback, (timer_proc == Qnil ? NULL : s_Molecule_CallSubProcessAsync_TimerCallback), fpout, fperr);
+
+ if (rb_obj_is_kind_of(cmd, rb_cString)) {
+ cmdargv = calloc(sizeof(cmdargv[0]), 3);
+ cmdargv[0] = StringValuePtr(cmd);
+ cmdargv[1] = "";
+ cmdargv[2] = NULL;
+ } else {
+ cmd = rb_ary_to_ary(cmd);
+ cmdargv = calloc(sizeof(cmdargv[0]), RARRAY_LEN(cmd) + 1);
+ for (n = 0; n < RARRAY_LEN(cmd); n++) {
+ cmdargv[n] = StringValuePtr(RARRAY_PTR(cmd)[n]);
+ }
+ cmdargv[n] = NULL;
+ }
+ n = MoleculeCallback_callSubProcessAsync(mol, cmdargv, s_Molecule_CallSubProcessAsync_EndCallback, (timer_proc == Qnil ? NULL : s_Molecule_CallSubProcessAsync_TimerCallback), fpout, fperr);
if (fpout != NULL && fpout != (FILE *)1)
fclose(fpout);
if (fperr != NULL && fperr != (FILE *)1)
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_method(rb_cMolecule, "set_molecule", s_Molecule_SetMolecule, 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, "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, "get_graphic_point", s_Molecule_GetGraphicPoint, -1);
+ rb_define_method(rb_cMolecule, "set_graphic_point", s_Molecule_SetGraphicPoint, -1);
+ rb_define_alias(rb_cMolecule, "get_graphic_points", "get_graphic_point");
+ rb_define_alias(rb_cMolecule, "set_graphic_points", "set_graphic_point");
+ rb_define_method(rb_cMolecule, "get_graphic_color", s_Molecule_SetGraphicColor, 1);
rb_define_method(rb_cMolecule, "set_graphic_color", s_Molecule_SetGraphicColor, 2);
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, "clear_surface", s_Molecule_ClearSurface, 0);
+ rb_define_method(rb_cMolecule, "show_surface", s_Molecule_ShowSurface, 0);
+ rb_define_method(rb_cMolecule, "hide_surface", s_Molecule_HideSurface, 0);
rb_define_method(rb_cMolecule, "create_surface", s_Molecule_CreateSurface, -1);
rb_define_method(rb_cMolecule, "set_surface_attr", s_Molecule_SetSurfaceAttr, 1);
rb_define_method(rb_cMolecule, "nelpots", s_Molecule_NElpots, 0);
rb_define_method(rb_cMolecule, "elpot", s_Molecule_Elpot, 1);
rb_define_method(rb_cMolecule, "clear_basis_set", s_Molecule_ClearBasisSet, 0);
- rb_define_method(rb_cMolecule, "add_gaussian_orbital_shell", s_Molecule_AddGaussianOrbitalShell, 3);
+ rb_define_method(rb_cMolecule, "add_gaussian_orbital_shell", s_Molecule_AddGaussianOrbitalShell, -1);
rb_define_method(rb_cMolecule, "add_gaussian_primitive_coefficients", s_Molecule_AddGaussianPrimitiveCoefficients, 3);
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);
rb_define_alias(rb_cAtomRef, "set_attr", "[]=");
rb_define_method(rb_cAtomRef, "[]", s_AtomRef_GetAttr, 1);
rb_define_alias(rb_cAtomRef, "get_attr", "[]");
- s_SetAtomAttrString = rb_str_new2("set_atom_attr");
+ s_SetAtomAttrString = Ruby_NewEncodedStringValue2("set_atom_attr");
rb_global_variable(&s_SetAtomAttrString);
rb_define_method(rb_cAtomRef, "molecule", s_AtomRef_GetMolecule, 0);
rb_define_method(rb_cAtomRef, "==", s_AtomRef_Equal, 1);
rb_define_method(rb_mKernel, "play_sound", s_Kernel_PlaySound, -1);
rb_define_method(rb_mKernel, "stop_sound", s_Kernel_StopSound, 0);
rb_define_method(rb_mKernel, "export_to_clipboard", s_Kernel_ExportToClipboard, 1);
+ rb_define_method(rb_mKernel, "hartree_to_kcal", s_Kernel_HartreeToKcal, 1);
+ rb_define_method(rb_mKernel, "hartree_to_kj", s_Kernel_HartreeToKJ, 1);
+ rb_define_method(rb_mKernel, "kcal_to_hartree", s_Kernel_KcalToHartree, 1);
+ rb_define_method(rb_mKernel, "kj_to_hartree", s_Kernel_KJToHartree, 1);
+ rb_define_method(rb_mKernel, "bohr_to_angstrom", s_Kernel_BohrToAngstrom, 1);
+ rb_define_method(rb_mKernel, "angstrom_to_bohr", s_Kernel_AngstromToBohr, 1);
/* class IO */
rb_define_method(rb_cIO, "gets_any_eol", s_IO_gets_any_eol, 0);
g_RubyID_call = rb_intern("call");
s_InitMOInfoKeys();
+
+ /* Symbols for graphics */
+ s_LineSym = ID2SYM(rb_intern("line"));
+ s_PolySym = ID2SYM(rb_intern("poly"));
+ s_CylinderSym = ID2SYM(rb_intern("cylinder"));
+ s_ConeSym = ID2SYM(rb_intern("cone"));
+ s_EllipsoidSym = ID2SYM(rb_intern("ellipsoid"));
+}
+
+#pragma mark ====== Interface with RubyDialog class ======
+
+RubyValue
+RubyDialogCallback_parentModule(void)
+{
+ return (RubyValue)rb_mMolby;
}
#pragma mark ====== External functions ======
#else
asprintf(&scr, "#coding:utf-8\n%s", (char *)ptr[0]);
#endif
- sval = rb_str_new2(scr);
+ sval = Ruby_NewEncodedStringValue2(scr);
free(scr);
- fnval = rb_str_new2("(eval)");
+ fnval = Ruby_NewEncodedStringValue2("(eval)");
lnval = INT2FIX(0);
} else {
- sval = rb_str_new2((char *)ptr[0]);
+ sval = Ruby_NewEncodedStringValue2((char *)ptr[0]);
fnval = Ruby_NewFileStringValue((char *)ptr[2]);
lnval = INT2FIX(1);
}
return retval;
}
+/* For debug */
+char *
+Ruby_inspectValue(RubyValue value)
+{
+ int status;
+ static char buf[256];
+ VALUE val = (VALUE)value;
+ gMolbyRunLevel++;
+ val = rb_protect(rb_inspect, val, &status);
+ gMolbyRunLevel--;
+ if (status == 0) {
+ char *str = StringValuePtr(val);
+ strncpy(buf, str, sizeof(buf) - 1);
+ buf[sizeof(buf) - 1] = 0;
+ } else {
+ snprintf(buf, sizeof(buf), "Error status = %d", status);
+ }
+ return buf;
+}
+
int
-Molby_showRubyValue(RubyValue value, char **outValueString)
+Ruby_showValue(RubyValue value, char **outValueString)
{
VALUE val = (VALUE)value;
if (gMolbyIsCheckingInterrupt) {
}
void
-Molby_showError(int status)
+Ruby_showError(int status)
{
static const int tag_raise = 6;
+ char *main_message = "Molby script error";
char *msg = NULL, *msg2;
VALUE val, backtrace;
int interrupted = 0;
+ int exit_status = -1;
if (status == tag_raise) {
VALUE errinfo = rb_errinfo();
VALUE eclass = CLASS_OF(errinfo);
if (eclass == rb_eInterrupt) {
- msg = "Interrupt";
+ main_message = "Molby script interrupted";
+ msg = "Interrupt";
interrupted = 1;
- }
+ } else if (eclass == rb_eSystemExit) {
+ main_message = "Molby script exit";
+ interrupted = 2;
+ val = rb_eval_string_protect("$!.status", &status);
+ if (status == 0) {
+ exit_status = NUM2INT(rb_Integer(val));
+ asprintf(&msg, "Molby script exit with status %d", exit_status);
+ } else {
+ asprintf(&msg, "Molby script exit with unknown status");
+ }
+ }
}
gMolbyRunLevel++;
- backtrace = rb_eval_string_protect("$backtrace = $!.backtrace.join(\"\\n\")", &status);
- if (msg == NULL) {
- val = rb_eval_string_protect("$!.to_s", &status);
- if (status == 0)
- msg = RSTRING_PTR(val);
- else msg = "(message not available)";
- }
- asprintf(&msg2, "%s\n%s", msg, RSTRING_PTR(backtrace));
- MyAppCallback_messageBox(msg2, (interrupted == 0 ? "Molby script error" : "Molby script interrupted"), 0, 3);
+ if (exit_status != 0) {
+ backtrace = rb_eval_string_protect("$backtrace = $!.backtrace.join(\"\\n\")", &status);
+ if (msg == NULL) {
+ val = rb_eval_string_protect("$!.to_s", &status);
+ if (status == 0)
+ msg = RSTRING_PTR(val);
+ else
+ msg = "(message not available)";
+ }
+ asprintf(&msg2, "%s\n%s", msg, RSTRING_PTR(backtrace));
+ } else {
+ msg2 = strdup(msg);
+ }
+ MyAppCallback_messageBox(msg2, main_message, 0, 3);
free(msg2);
+ if (interrupted == 2) {
+ free(msg);
+ if (!gUseGUI && exit_status == 0)
+ exit(0); // Capture exit(0) here and force exit
+ }
gMolbyRunLevel--;
}
-char *
-Molby_getDescription(void)
+/* Wrapper function for rb_load_protect or rb_eval_string_protect. Used only in non-GUI mode. */
+int
+Molby_loadScript(const char *script, int from_file)
+{
+ int status;
+ gMolbyRunLevel++;
+ if (from_file)
+ rb_load_protect(Ruby_NewEncodedStringValue2(script), 0, &status);
+ else
+ rb_eval_string_protect(script, &status);
+ gMolbyRunLevel--;
+ return status;
+}
+
+void
+Molby_getDescription(char **versionString, char **auxString)
{
extern const char *gVersionString, *gCopyrightString;
extern int gRevisionNumber;
extern char *gLastBuildString;
- char *s;
+ char *s1, *s2;
char *revisionString;
if (gRevisionNumber > 0) {
asprintf(&revisionString, ", revision %d", gRevisionNumber);
} else revisionString = "";
- asprintf(&s,
- "Molby %s%s\n%s\nLast compile: %s\n"
-#if !defined(__CMDMAC__)
- "\nIncluding:\n"
- "%s"
+
+ asprintf(&s1, "%s %s%s\n%s\nLast compile: %s\n",
+#if defined(__WXMSW__)
+ #if TARGET_ARCH == 64
+ "Molby (64 bit)",
+ #else
+ "Molby (32 bit)",
+ #endif
#else
- "Including "
-#endif
- "ruby %s, http://www.ruby-lang.org/\n"
- "%s\n"
- "FFTW 3.3.2, http://www.fftw.org/\n"
- " Copyright (C) 2003, 2007-11 Matteo Frigo\n"
- " and Massachusetts Institute of Technology",
- gVersionString, revisionString, gCopyrightString, gLastBuildString,
-#if !defined(__CMDMAC__)
- MyAppCallback_getGUIDescriptionString(),
+ "Molby",
#endif
- gRubyVersion, gRubyCopyright);
+ gVersionString, revisionString, gCopyrightString, gLastBuildString);
+ if (gUseGUI) {
+ asprintf(&s2,
+ "\nIncluding:\n"
+ "%s"
+ "ruby %s, http://www.ruby-lang.org/\n"
+ "%s\n"
+ "FFTW 3.3.2, http://www.fftw.org/\n"
+ " Copyright (C) 2003, 2007-11 Matteo Frigo"
+ " and Massachusetts Institute of Technology",
+ MyAppCallback_getGUIDescriptionString(),
+ gRubyVersion, gRubyCopyright);
+ } else {
+ asprintf(&s2,
+ "Including "
+ "ruby %s, http://www.ruby-lang.org/\n"
+ "%s\n"
+ "FFTW 3.3.2, http://www.fftw.org/\n"
+ " Copyright (C) 2003, 2007-11 Matteo Frigo"
+ " and Massachusetts Institute of Technology",
+ gRubyVersion, gRubyCopyright);
+
+ }
if (revisionString[0] != 0)
free(revisionString);
- return s;
+ if (versionString != NULL)
+ *versionString = s1;
+ if (auxString != NULL)
+ *auxString = s2;
}
void
{
gRubyVersion = strdup(ruby_version);
asprintf(&gRubyCopyright, "%sCopyright (C) %d-%d %s",
-#if defined(__CMDMAC__)
- "",
-#else
" ", /* Indent for displaying in About dialog */
-#endif
RUBY_BIRTH_YEAR, RUBY_RELEASE_YEAR, RUBY_AUTHOR);
}
}
} */
-#if defined(__CMDMAC__)
- wbuf = Molby_getDescription();
- printf("%s\n", wbuf);
- free(wbuf);
-#endif
+ if (!gUseGUI) {
+ char *wbuf2;
+ Molby_getDescription(&wbuf, &wbuf2);
+ MyAppCallback_showScriptMessage("%s\n%s\n", wbuf, wbuf2);
+ free(wbuf);
+ free(wbuf2);
+ }
/* Read atom display parameters */
if (ElementParameterInitialize("element.par", &wbuf) != 0) {
-#if defined(__CMDMAC__)
- fprintf(stderr, "%s\n", wbuf);
-#else
- MyAppCallback_setConsoleColor(1);
- MyAppCallback_showScriptMessage("%s", wbuf);
- MyAppCallback_setConsoleColor(0);
-#endif
+ MyAppCallback_setConsoleColor(1);
+ MyAppCallback_showScriptMessage("%s", wbuf);
+ MyAppCallback_setConsoleColor(0);
free(wbuf);
}
/* Read default parameters */
ParameterReadFromFile(gBuiltinParameters, "default.par", &wbuf, NULL);
if (wbuf != NULL) {
-#if defined(__CMDMAC__)
- fprintf(stderr, "%s\n", wbuf);
-#else
- MyAppCallback_setConsoleColor(1);
- MyAppCallback_showScriptMessage("%s", wbuf);
- MyAppCallback_setConsoleColor(0);
-#endif
+ MyAppCallback_setConsoleColor(1);
+ MyAppCallback_showScriptMessage("%s", wbuf);
+ MyAppCallback_setConsoleColor(0);
free(wbuf);
}
/* Initialize Ruby interpreter */
#if __WXMSW__
- {
+ if (gUseGUI) {
/* On Windows, fileno(stdin|stdout|stderr) returns -2 and
it causes rb_bug() (= fatal error) during ruby_init().
As a workaround, these standard streams are reopend as
ruby_init();
{
- extern void Init_shift_jis(void), Init_trans_japanese_sjis(void);
- Init_shift_jis();
- Init_trans_japanese_sjis();
+ /* Initialize CP932/Windows-31J encodings */
+ extern void Init_shift_jis(void), Init_windows_31j(void), Init_trans_japanese_sjis(void);
+ extern int rb_enc_alias(const char *, const char *);
+ Init_shift_jis();
+ Init_windows_31j();
+ Init_trans_japanese_sjis();
+ rb_enc_alias("CP932", "Windows-31J");
+ }
+
+#if defined(__WXMSW__)
+ {
+ /* Set default external encoding */
+ /* The following snippet is taken from encoding.c */
+ extern void rb_enc_set_default_external(VALUE encoding);
+ char cp[sizeof(int) * 8 / 3 + 22];
+ int status;
+ VALUE enc;
+ snprintf(cp, sizeof cp, "Encoding.find('CP%d')", AreFileApisANSI() ? GetACP() : GetOEMCP());
+ enc = rb_eval_string_protect(cp, &status);
+ if (status == 0 && !NIL_P(enc)) {
+ rb_enc_set_default_external(enc);
+ }
}
-
+#endif
+
/* Initialize loadpath; the specified directory, "lib" subdirectory, and "." */
ruby_incpush(".");
asprintf(&libpath, "%s%clib", dir, PATH_SEPARATOR);
/* Define Molby classes */
Init_Molby();
- RubyDialogInitClass();
+ if (gUseGUI)
+ RubyDialogInitClass();
rb_define_const(rb_mMolby, "ResourcePath", val);
val = Ruby_NewFileStringValue(dir);
rb_define_const(rb_mMolby, "DocumentDirectory", val);
free(p);
-#if defined(__CMDMAC__)
- rb_define_const(rb_mMolby, "HasGUI", Qfalse);
-#else
- rb_define_const(rb_mMolby, "HasGUI", Qtrue);
-#endif
-
-#if !__CMDMAC__
-
- /* Create objects for stdout and stderr */
- val = rb_funcall(rb_cObject, rb_intern("new"), 0);
- rb_define_singleton_method(val, "write", s_StandardOutput, 1);
- rb_define_singleton_method(val, "flush", s_FlushConsoleOutput, 0);
- rb_gv_set("$stdout", val);
- val = rb_funcall(rb_cObject, rb_intern("new"), 0);
- rb_define_singleton_method(val, "write", s_StandardErrorOutput, 1);
- rb_define_singleton_method(val, "flush", s_FlushConsoleOutput, 0);
- rb_gv_set("$stderr", val);
-
- /* Create objects for stdin */
- val = rb_funcall(rb_cObject, rb_intern("new"), 0);
- rb_define_singleton_method(val, "gets", s_StandardInputGets, -1);
- rb_define_singleton_method(val, "readline", s_StandardInputGets, -1);
- rb_define_singleton_method(val, "method_missing", s_StandardInputMethodMissing, -1);
- rb_gv_set("$stdin", val);
-
-#endif
+ if (gUseGUI)
+ rb_define_const(rb_mMolby, "HasGUI", Qtrue);
+ else
+ rb_define_const(rb_mMolby, "HasGUI", Qfalse);
+
+ {
+ /* Create objects for stdout and stderr */
+ val = rb_funcall(rb_cObject, rb_intern("new"), 0);
+ rb_define_singleton_method(val, "write", s_StandardOutput, 1);
+ rb_define_singleton_method(val, "flush", s_FlushConsoleOutput, 0);
+ rb_gv_set("$stdout", val);
+ val = rb_funcall(rb_cObject, rb_intern("new"), 0);
+ rb_define_singleton_method(val, "write", s_StandardErrorOutput, 1);
+ rb_define_singleton_method(val, "flush", s_FlushConsoleOutput, 0);
+ rb_gv_set("$stderr", val);
+
+ /* Create objects for stdin */
+ val = rb_funcall(rb_cObject, rb_intern("new"), 0);
+ rb_define_singleton_method(val, "gets", s_StandardInputGets, -1);
+ rb_define_singleton_method(val, "readline", s_StandardInputGets, -1);
+ rb_define_singleton_method(val, "method_missing", s_StandardInputMethodMissing, -1);
+ rb_gv_set("$stdin", val);
+ }
/* Global variable to hold error information */
rb_define_variable("$backtrace", &gMolbyBacktrace);
gScriptMenuCommands = rb_ary_new();
gScriptMenuEnablers = rb_ary_new();
-#if !__CMDMAC__
- /* Register interrupt check code */
- rb_add_event_hook(s_Event_Callback, RUBY_EVENT_ALL, Qnil);
-#endif
-
-#if !__CMDMAC__
- /* Start interval timer (for periodic polling of interrupt); firing every 50 msec */
- s_SetIntervalTimer(0, 50);
-#endif
+ if (gUseGUI) {
+ /* Register interrupt check code */
+ rb_add_event_hook(s_Event_Callback, RUBY_EVENT_ALL, Qnil);
+ /* Start interval timer (for periodic polling of interrupt); firing every 50 msec */
+ s_SetIntervalTimer(0, 50);
+ }
/* Read the startup script */
if (script != NULL && script[0] != 0) {
MyAppCallback_showScriptMessage("Evaluating %s...\n", script);
gMolbyRunLevel++;
- rb_load_protect(rb_str_new2(script), 0, &status);
+ rb_load_protect(Ruby_NewEncodedStringValue2(script), 0, &status);
gMolbyRunLevel--;
if (status != 0)
- Molby_showError(status);
+ Ruby_showError(status);
else
MyAppCallback_showScriptMessage("Done.\n");
}
OSDN Git Service
