const char *gMolActionSetAtomForces = "atomForce:GV";
const char *gMolActionInsertFrames = "insertFrames:GVV";
const char *gMolActionRemoveFrames = "removeFrames:G";
+const char *gMolActionReorderFrames = "reorderFrames:I";
const char *gMolActionSetProperty = "setProperty:iGD";
const char *gMolActionSetSelection = "selection:G";
const char *gMolActionChangeResidueNumber = "changeResSeq:Gi";
}
static int
+s_MolActionReorderFrames(Molecule *mol, MolAction *action, MolAction **actp)
+{
+ Int *ip, *ip2, i, n;
+
+ ip = action->args[0].u.arval.ptr;
+
+ /* Prepare undo data (and error check) */
+ n = MoleculeGetNumberOfFrames(mol);
+ ip2 = (Int *)malloc(sizeof(Int) * n);
+ for (i = 0; i < n; i++) {
+ if (ip[i] < 0 || ip[i] >= n) {
+ free(ip2);
+ return -1;
+ }
+ ip2[ip[i]] = i;
+ }
+
+ MoleculeReorderFrames(mol, ip);
+
+ *actp = MolActionNew(gMolActionReorderFrames, n, ip2);
+ (*actp)->frame = mol->cframe;
+ free(ip2);
+
+ return 0;
+}
+
+static int
s_MolActionSetProperty(Molecule *mol, MolAction *action, MolAction **actp)
{
Int idx, n1, no_undo;
if ((result = s_MolActionRemoveFrames(mol, action, &act2)) != 0)
return result;
needsSymmetryAmendment = 1;
+ } else if (strcmp(action->name, gMolActionReorderFrames) == 0) {
+ if ((result = s_MolActionReorderFrames(mol, action, &act2)) != 0)
+ return result;
} else if (strcmp(action->name, gMolActionSetProperty) == 0) {
if ((result = s_MolActionSetProperty(mol, action, &act2)) != 0)
return result;
extern const char *gMolActionSetAtomForces;
extern const char *gMolActionInsertFrames;
extern const char *gMolActionRemoveFrames;
+extern const char *gMolActionReorderFrames;
extern const char *gMolActionSetProperty;
extern const char *gMolActionSetSelection;
extern const char *gMolActionChangeResidueNumber;
ParameterRelease(mp->par);
mp->par = NULL;
}
- if (mp->bset != NULL) {
- BasisSetRelease(mp->bset);
- mp->bset = NULL;
- }
if (mp->atoms != NULL) {
for (i = 0; i < mp->natoms; i++)
AtomClean(mp->atoms + i);
mp->bset = NULL;
}
if (mp->mcube != NULL) {
- free(mp->mcube->dp);
- free(mp->mcube->radii);
- free(mp->mcube->c[0].cubepoints);
- free(mp->mcube->c[0].triangles);
- free(mp->mcube->c[1].cubepoints);
- free(mp->mcube->c[1].triangles);
- free(mp->mcube);
+ MoleculeDeallocateMCube(mp->mcube);
mp->mcube = NULL;
}
if (mp->molprops != NULL) {
if (fn > 1) {
for (i = 0; i < mp->natoms; i++) {
ap = ATOM_AT_INDEX(mp->atoms, i);
- ap->frames = (Vector *)malloc(sizeof(Vector) * fn);
+ NewArray(&ap->frames, &ap->nframes, sizeof(Vector), fn);
if (ap->frames == NULL)
goto panic;
- ap->nframes = fn;
for (j = 0; j < fn; j++)
ap->frames[j] = frames[mp->natoms * j + i];
}
return 0;
}
-/* Allocate BasisSet record. rflag: UHF, 0; RHF, 1; ROHF, 2
+/* Allocate BasisSet record. rflag: 0, UHF; 1, RHF; 2, ROHF; -1, clear
ne_alpha: number of alpha electrons, ne_beta: number of beta electrons
The natoms and pos are copied from mol. */
int
MoleculeAllocateBasisSetRecord(Molecule *mol, Int rflag, Int ne_alpha, Int ne_beta)
{
BasisSet *bset;
-/* int i;
- Atom *ap; */
if (mol == NULL || mol->natoms == 0)
return -1; /* Molecule is empty */
+ if (rflag < 0) {
+ if (mol->bset != NULL) {
+ BasisSetRelease(mol->bset);
+ mol->bset = NULL;
+ }
+ return 0;
+ }
bset = mol->bset;
if (bset == NULL) {
bset = mol->bset = (BasisSet *)calloc(sizeof(BasisSet), 1);
for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
if (ap->nframes == 0)
continue;
- ap->frames = (Vector *)malloc(sizeof(Vector) * ap->nframes);
+ n = ap->nframes;
+ ap->frames = NULL;
+ ap->nframes = 0;
+ NewArray(&ap->frames, &ap->nframes, sizeof(Vector), n);
if (ap->frames == NULL)
goto out_of_memory;
memmove(ap->frames, ptr, sizeof(Vector) * ap->nframes);
/* Expand ap->frames for all atoms */
for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
- if (ap->frames == NULL)
- vp = (Vector *)calloc(sizeof(Vector), n_new);
- else
- vp = (Vector *)realloc(ap->frames, sizeof(Vector) * n_new);
- if (vp == NULL) {
+ Int n = ap->nframes;
+ AssignArray(&ap->frames, &ap->nframes, sizeof(Vector), n_new - 1, NULL);
+ if (ap->frames == NULL) {
__MoleculeUnlock(mp);
return -1;
}
- for (j = ap->nframes; j < n_new; j++)
- vp[j] = ap->r;
- ap->frames = vp;
+ for (j = n; j < n_new; j++)
+ ap->frames[j] = ap->r;
}
if (mp->cell != NULL) {
j = mp->nframe_cells;
mult = sizeof(Vector);
vp = ap->frames;
if (vp == NULL) {
- ap->frames = vp = (Vector *)calloc(sizeof(Vector), n_old);
- if (vp == NULL) {
+ NewArray(&ap->frames, &ap->nframes, sizeof(Vector), n_old);
+ if (ap->frames == NULL) {
__MoleculeUnlock(mp);
return -1;
}
+ vp = ap->frames;
}
old_count = ap->nframes;
} else {
prp->propvals = (Double *)realloc(prp->propvals, sizeof(Double));
}
} else {
- if (i < mp->natoms)
- ap->frames = (Vector *)realloc(ap->frames, sizeof(Vector) * s);
- else if (i == mp->natoms) {
+ if (i < mp->natoms) {
+ AssignArray(&ap->frames, &ap->nframes, sizeof(Vector), s - 1, NULL);
+ ap->nframes = s;
+ } else if (i == mp->natoms) {
AssignArray(&mp->frame_cells, &mp->nframe_cells, sizeof(Vector) * 4, s - 1, NULL);
mp->nframe_cells = s;
} else {
/* Write the current coordinate back to the frame array */
ap->frames[cframe] = ap->r;
}
- if (frame != cframe && frame >= 0 && frame < ap->nframes) {
+ if ((frame != cframe || copyback == 0) && frame >= 0 && frame < ap->nframes) {
/* Read the coordinate from the frame array */
ap->r = ap->frames[frame];
modified = 1;
vp[3] = mp->cell->origin;
}
/* Set the cell from the frame array */
- if (frame != cframe && frame >= 0 && frame < mp->nframe_cells) {
+ if ((frame != cframe || copyback == 0) && frame >= 0 && frame < mp->nframe_cells) {
MoleculeSetPeriodicBox(mp, &mp->frame_cells[frame * 4], &mp->frame_cells[frame * 4 + 1], &mp->frame_cells[frame * 4 + 2], &mp->frame_cells[frame * 4 + 3], mp->cell->flags, 0);
modified = 1;
MoleculeAmendBySymmetry(mp, NULL, NULL, NULL);
return nframes;
}
+int
+MoleculeReorderFrames(Molecule *mp, const Int *old_idx)
+{
+ Int *ip, i, j, n, nframes;
+ Double *dp;
+ Atom *ap;
+ MolProp *prp;
+ if (mp == NULL || old_idx == NULL)
+ return 0;
+ nframes = MoleculeGetNumberOfFrames(mp);
+ MoleculeFlushFrames(mp);
+ ip = (Int *)malloc(sizeof(Int) * nframes);
+ if (ip == NULL)
+ return -1; /* Out of memory */
+ memset(ip, 0, sizeof(Int) * nframes);
+ /* Check the argument */
+ for (i = 0; i < nframes; i++) {
+ j = old_idx[i];
+ if (j < 0 || j >= nframes || ip[j] != 0) {
+ free(ip);
+ return -2; /* Bad argument */
+ }
+ ip[j] = 1;
+ }
+ free(ip);
+ dp = (Double *)malloc(sizeof(Double) * nframes * 12);
+ for (i = 0, ap = mp->atoms, prp = mp->molprops; i <= mp->natoms + mp->nmolprops; i++) {
+ for (j = 0; j < nframes; j++) {
+ n = old_idx[j];
+ if (i < mp->natoms) {
+ ((Vector *)dp)[j] = (n < ap->nframes ? ap->frames[n] : ap->r);
+ } else if (i == mp->natoms) {
+ if (mp->cell != NULL) {
+ if (n < mp->nframe_cells && mp->frame_cells != NULL)
+ memmove(dp + j * 12, mp->frame_cells + n * 4, sizeof(Vector) * 4);
+ else {
+ ((Vector *)dp)[j * 4] = mp->cell->axes[0];
+ ((Vector *)dp)[j * 4] = mp->cell->axes[1];
+ ((Vector *)dp)[j * 4] = mp->cell->axes[2];
+ ((Vector *)dp)[j * 4] = mp->cell->origin;
+ }
+ }
+ } else {
+ dp[j] = prp->propvals[n];
+ }
+ }
+ for (j = 0; j < nframes; j++) {
+ if (i < mp->natoms) {
+ if (ap->nframes <= j)
+ AssignArray(&ap->frames, &ap->nframes, sizeof(Vector), nframes - 1, NULL);
+ ap->frames[j] = ((Vector *)dp)[j];
+ } else if (i == mp->natoms) {
+ if (mp->cell != NULL) {
+ AssignArray(&mp->frame_cells, &mp->nframe_cells, sizeof(Vector) * 4, nframes - 1, NULL);
+ memmove(mp->frame_cells + j * 4, dp + j * 12, sizeof(Vector) * 4);
+ }
+ } else {
+ prp->propvals[j] = dp[j];
+ }
+ }
+ if (i < mp->natoms)
+ ap = ATOM_NEXT(ap);
+ else if (i > mp->natoms)
+ prp++;
+ }
+ free(dp);
+ MoleculeSelectFrame(mp, mp->cframe, 0);
+ return 0;
+}
+
#pragma mark ====== Molecule Propeties ======
int
return -10000;
AssignArray(&mp->molprops, &mp->nmolprops, sizeof(MolProp), mp->nmolprops, prp);
free(prp);
- return mp->nmolprops;
+ return mp->nmolprops - 1;
}
int
return retval;
}
+
+void
+MoleculeDeallocateMCube(MCube *mcube)
+{
+ free(mcube->dp);
+ free(mcube->radii);
+ free(mcube->c[0].cubepoints);
+ free(mcube->c[0].triangles);
+ free(mcube->c[1].cubepoints);
+ free(mcube->c[1].triangles);
+ free(mcube);
+}
int strlen_limit(const char *s, int limit);
+void BasisSetRelease(BasisSet *bset);
+
Molecule *MoleculeNew(void);
int MoleculeLoadFile(Molecule *mp, const char *fname, const char *ftype, char **errbuf);
int MoleculeLoadPsfFile(Molecule *mp, const char *fname, char **errbuf);
int MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame, Vector *outFrameCell);
int MoleculeSelectFrame(Molecule *mp, int frame, int copyback);
int MoleculeFlushFrames(Molecule *mp);
+int MoleculeReorderFrames(Molecule *mp, const Int *old_idx);
int MoleculeCreateProperty(Molecule *mp, const char *name);
int MoleculeLookUpProperty(Molecule *mp, const char *name);
MCube *MoleculeClearMCube(Molecule *mol, Int nx, Int ny, Int nz, const Vector *origin, Double dx, Double dy, Double dz);
int MoleculeUpdateMCube(Molecule *mol, int idn);
+void MoleculeDeallocateMCube(MCube *mcube);
extern char *gMoleculePasteboardType;
extern char *gParameterPasteboardType;
/*
* 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:
* get_coord_from_frame(index, group = nil)
*
* Copy the coordinates from the indicated frame. If group is specified, only the specified atoms
}
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:
/*
* call-seq:
+ * clear_basis_set
+ *
+ * Clear the existing basis set info. All gaussian coefficients, MO energies and coefficients,
+ * cube and marching cube information are discarded. This operation is _not_ undoable!
+ */
+static VALUE
+s_Molecule_ClearBasisSet(VALUE self)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol != NULL) {
+ if (mol->bset != NULL) {
+ BasisSetRelease(mol->bset);
+ mol->bset = NULL;
+ }
+ if (mol->mcube != NULL) {
+ MoleculeDeallocateMCube(mol->mcube);
+ mol->mcube = NULL;
+ }
+ }
+ return self;
+}
+
+/*
+ * call-seq:
* add_gaussian_orbital_shell(sym, nprims, atom_index)
*
* To be used internally. Add a gaussian orbital shell with symmetry code, number of primitives,
/*
* call-seq:
- * mo_type
+ * clear_mo_coefficients
*
- * Returns either "RHF", "UHF", or "ROHF". If no MO info is present, returns nil.
+ * Clear the existing MO coefficients.
*/
static VALUE
-s_Molecule_MOType(VALUE self)
+s_Molecule_ClearMOCoefficients(VALUE self)
{
Molecule *mol;
- Data_Get_Struct(self, Molecule, mol);
- if (mol != NULL && mol->bset != NULL) {
- const char *s;
- int rflag = mol->bset->rflag;
- if (rflag == 0)
- s = "UHF";
- else if (rflag == 2)
- s = "ROHF";
- else s = "RHF";
- return rb_str_new2(s);
- } else return Qnil;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->bset != NULL) {
+ if (mol->bset->moenergies != NULL) {
+ free(mol->bset->moenergies);
+ mol->bset->moenergies = NULL;
+ }
+ if (mol->bset->mo != NULL) {
+ free(mol->bset->mo);
+ mol->bset->mo = NULL;
+ }
+ mol->bset->nmos = 0;
+ }
+ return self;
}
/*
return rb_float_new(energy);
}
+static VALUE sTypeSym, sAlphaSym, sBetaSym;
+
+static inline void
+s_InitMOInfoKeys(void)
+{
+ if (sTypeSym == 0) {
+ sTypeSym = ID2SYM(rb_intern("type"));
+ sAlphaSym = ID2SYM(rb_intern("alpha"));
+ sBetaSym = ID2SYM(rb_intern("beta"));
+ }
+}
+
+/*
+ * call-seq:
+ * set_mo_info(hash)
+ *
+ * Set the MO info. hash keys: :type=>"RHF"|"UHF"|"ROHF",
+ * :alpha=>integer, :beta=>integer
+ */
+static VALUE
+s_Molecule_SetMOInfo(VALUE self, VALUE hval)
+{
+ Molecule *mol;
+ VALUE aval;
+ Int rflag, na, nb, n;
+ char *s;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->bset != NULL) {
+ rflag = mol->bset->rflag;
+ na = mol->bset->ne_alpha;
+ nb = mol->bset->ne_beta;
+ } else {
+ rflag = 1;
+ na = 0;
+ nb = 0;
+ }
+ if (hval != Qnil) {
+ if ((aval = rb_hash_aref(hval, sTypeSym)) != Qnil) {
+ s = StringValuePtr(aval);
+ if (strcasecmp(s, "RHF") == 0)
+ rflag = 1;
+ else if (strcasecmp(s, "UHF") == 0)
+ rflag = 0;
+ else if (strcasecmp(s, "ROHF") == 0)
+ rflag = 2;
+ }
+ if ((aval = rb_hash_aref(hval, sAlphaSym)) != Qnil) {
+ n = NUM2INT(rb_Integer(aval));
+ if (n >= 0)
+ na = n;
+ }
+ if ((aval = rb_hash_aref(hval, sBetaSym)) != Qnil) {
+ n = NUM2INT(rb_Integer(aval));
+ if (n >= 0)
+ nb = n;
+ }
+ MoleculeAllocateBasisSetRecord(mol, rflag, na, nb);
+ }
+ return self;
+}
+
+/*
+ * call-seq:
+ * get_mo_info(key)
+ *
+ * Get the MO info. The key is as described in set_mo_info.
+ */
+static VALUE
+s_Molecule_GetMOInfo(VALUE self, VALUE kval)
+{
+ Molecule *mol;
+ Data_Get_Struct(self, Molecule, mol);
+ if (mol->bset == NULL)
+ rb_raise(rb_eMolbyError, "No MO information is defined");
+ 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");
+ default: return rb_str_to_str(INT2NUM(mol->bset->rflag));
+ }
+ } else if (kval == sAlphaSym) {
+ return INT2NUM(mol->bset->ne_alpha);
+ } else if (kval == sBetaSym) {
+ return INT2NUM(mol->bset->ne_beta);
+ } else return Qnil;
+}
+
+/*
+ * call-seq:
+ * mo_type
+ *
+ * Returns either "RHF", "UHF", or "ROHF". If no MO info is present, raises exception.
+ */
+static VALUE
+s_Molecule_MOType(VALUE self)
+{
+ return s_Molecule_GetMOInfo(self, sTypeSym);
+}
+
+
+#if 0
/*
* call-seq:
* allocate_basis_set_record(rflag, ne_alpha, ne_beta)
rb_raise(rb_eMolbyError, "Unknown error");
return self;
}
+#endif
/*
* call-seq:
{
Molecule *mol;
VALUE rval, nval;
- int i, n;
+ int i;
Data_Get_Struct(self, Molecule, mol);
rval = rb_ary_new();
for (i = mol->nmolprops - 1; i >= 0; i--) {
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, "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_primitive_coefficients", s_Molecule_AddGaussianPrimitiveCoefficients, 3);
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, "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, "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, "search_equivalent_atoms", s_Molecule_SearchEquivalentAtoms, -1);
rb_define_method(rb_cMolecule, "create_pi_anchor", s_Molecule_CreatePiAnchor, -1);
s_ID_equal = rb_intern("==");
g_RubyID_call = rb_intern("call");
+
+ s_InitMOInfoKeys();
}
#pragma mark ====== External functions ======
end
end
+ def cmd_reverse_frames
+ n = nframes
+ return if n == 0
+ hash = Dialog.run("Reverse Frames") {
+ layout(2,
+ item(:text, :title=>"Start"),
+ item(:textfield, :width=>120, :tag=>"start", :value=>"0"),
+ item(:text, :title=>"End"),
+ item(:textfield, :width=>120, :tag=>"end", :value=>(n - 1).to_s))
+ }
+ if hash[:status] == 0
+ sframe = Integer(hash["start"])
+ eframe = Integer(hash["end"])
+ return if sframe > eframe
+ eframe = n - 1 if eframe >= n
+ frames = (0...sframe).to_a + (sframe..eframe).to_a.reverse + ((eframe + 1)...n).to_a
+ reorder_frames(frames)
+ end
+ end
+
def cmd_extra_properties
mol = self
- on_document_modified = lambda { |m|
- }
get_count = lambda { |it| mol.nframes }
- get_value = lambda { |it, row, col| sprintf("%.8g", mol.get_property(col, row)) rescue "" }
+ get_value = lambda { |it, row, col|
+ if col == 0
+ row.to_s
+ else
+ sprintf("%.8f", mol.get_property(col - 1, row)) rescue ""
+ end
+ }
Dialog.new("Extra Props:" + mol.name, nil, nil, :resizable=>true, :has_close_box=>true) {
- columns = mol.property_names.map { |name| [name, 80] }
+ names = nil
+ on_document_modified = lambda { |m|
+ if (n = m.property_names) != names
+ names = n
+ col = [["frame", 40]] + names.map { |nn| [nn, 120] }
+ item_with_tag("table")[:columns] = col
+ end
+ item_with_tag("table")[:refresh] = true
+ }
layout(1,
- item(:table, :width=>240, :height=>380, :flex=>[0,0,0,0,1,1], :tag=>"table",
- :columns=>columns,
+ item(:table, :width=>320, :height=>300, :flex=>[0,0,0,0,1,1], :tag=>"table",
:on_count=>get_count,
:on_get_value=>get_value),
:flex=>[0,0,0,0,1,1]
)
- set_min_size(480, 200)
+ set_min_size(320, 200)
listen(mol, "documentModified", on_document_modified)
listen(mol, "documentWillClose", lambda { |m| close } )
+ on_document_modified.call(mol)
show
}
end
register_menu("Sort by residue", :cmd_sort_by_residue, :non_empty)
register_menu("", "")
register_menu("Delete Frames...", :cmd_delete_frames, lambda { |m| m && m.nframes > 1 } )
+register_menu("Reverse Frames...", :cmd_reverse_frames, lambda { |m| m && m.nframes > 1 } )
register_menu("", "")
register_menu("Open Extra Properties Window...", :cmd_extra_properties, lambda { |m| m && m.property_names.count > 0 } )
#register_menu("cmd test", :cmd_test)
size = 0
lines = []
last_line = ""
-
+ search_mode = 0
+
# Callback procs
term_callback = lambda { |m, n|
msg = "GAMESS execution on #{inpbase} "
if line =~ /GEOMETRY SEARCH POINT NSERCH= *(\d+)/
nserch = $1.to_i
last_i = i
- elsif line =~ /NSERCH:/
+ search_mode = 1
+ energy = nil
+ elsif line =~ /CONSTRAINED OPTIMIZATION POINT/
+ search_mode = 2
+ energy = nil
+ elsif line =~ /POINT *([0-9]+) *ON THE REACTION PATH/
+ nserch = $1.to_i
+ last_i = i
+ elsif line =~ /NSERCH: *([0-9]+)/
# print line
if mol
- dummy, n, grad = line.match(/NSERCH:[^0-9]*([0-9]+).*GRAD[^0-9]*([-.0-9]+)/).to_a
+ n = $1
+ line =~ /E= +([-.0-9]+)/
+ energy = $1.to_f
+ line =~ /GRAD\. MAX= +([-.0-9]+)/
+ grad = $1
mol.show_text("Search: #{n}\nGradient: #{grad}")
+ mol.set_property("energy", energy)
end
last_i = i
+ elsif line =~ /TOTAL ENERGY += *([-.0-9]+)/
+ energy = $1
+ if mol && search_mode == 2
+ mol.show_text("Point: #{nserch}\nEnergy: #{energy}")
+ end
+ energy = energy.to_f
+ elsif line =~ /FINAL .* ENERGY IS +([-.0-9]+) *AFTER/
+ energy = $1.to_f
+ if mol && search_mode == 0
+ mol.set_property("energy", energy)
+ end
elsif nserch > 0 && line =~ /ddikick\.x/
last_i = -1
break
coords.push(Vector3D[x.to_f, y.to_f, z.to_f])
}
mol.create_frame([coords])
+ if search_mode == 1 && energy
+ mol.set_property("energy", energy)
+ end
mol.display
last_i = i + mol.natoms + 2
i = last_i # Skip the processed lines
alias :loadmdcrd :loadcrd
alias :savemdcrd :savecrd
+ def sub_load_gamess_log_basis_set(lines, lineno)
+ ln = 0
+ while (line = lines[ln])
+ ln += 1
+ break if line =~ /SHELL\s+TYPE\s+PRIMITIVE/
+ end
+ ln += 1
+ i = -1
+ nprims = 0
+ sym = -10 # undefined
+ ncomps = 0
+ while (line = lines[ln])
+ ln += 1
+ break if line =~ /TOTAL NUMBER OF BASIS SET/
+ if line =~ /^\s*$/
+ # End of one shell
+ add_gaussian_orbital_shell(sym, nprims, i)
+ nprims = 0
+ sym = -10
+ next
+ end
+ a = line.split
+ if a.length == 1
+ i += 1
+ ln += 1 # Skip the blank line
+ next
+ elsif a.length == 5 || a.length == 6
+ if sym == -10
+ case a[1]
+ when "S"
+ sym = 0; n = 1
+ when "P"
+ sym = 1; n = 3
+ when "L"
+ sym = -1; n = 4
+ when "D"
+ sym = 2; n = 6
+ when "F"
+ sym = 3; n = 10
+ when "G"
+ sym = 4; n = 15
+ else
+ raise MolbyError, "Unknown gaussian shell type at line #{lineno + ln}"
+ end
+ ncomps += n
+ end
+ if (a.length == 5 && sym == -1) || (a.length == 6 && sym != -1)
+ raise MolbyError, "Wrong format in gaussian shell information at line #{lineno + ln}"
+ end
+ exp = Float(a[3])
+ c = Float(a[4])
+ csp = Float(a[5] || 0.0)
+ add_gaussian_primitive_coefficients(exp, c, csp)
+ nprims += 1
+ else
+ raise MolbyError, "Error in reading basis set information at line #{lineno + ln}"
+ end
+ end
+ return ncomps
+ end
+
def sub_load_gamess_log(fp)
if natoms == 0
self.update_enabled = false
mes = "Loading GAMESS log file"
show_progress_panel(mes)
- energy = 0.0
+ energy = nil
ne_alpha = ne_beta = 0 # number of electrons
rflag = nil # 0, UHF; 1, RHF; 2, ROHF
mo_count = 0
+ search_mode = 0 # 0, no search; 1, optimize; 2, irc
ncomps = 0 # Number of AO terms per one MO (sum of the number of components over all shells)
alpha_beta = nil # Flag to read alpha/beta MO appropriately
+ nsearch = 0 # Search number for optimization
begin
- if nframes > 0
- create_frame
- frame = nframes - 1
- end
+ # if nframes > 0
+ # create_frame
+ # frame = nframes - 1
+ # end
n = 0
while 1
line = fp.gets
line.chomp!
if line =~ /ATOM\s+ATOMIC\s+COORDINATES/ || line =~ /COORDINATES OF ALL ATOMS ARE/
set_progress_message(mes + "\nReading atomic coordinates...")
- first_line = (line =~ /ATOMIC/)
+ if line =~ /ATOMIC/
+ first_line = true
+ if !new_unit
+ next # Skip initial atomic coordinates unless loading into an empty molecule
+ end
+ else
+ first_line = false
+ nsearch += 1
+ end
line = fp.gets # Skip one line
n = 0
coords = []
# }
# }
new_unit = false
- create_frame
+ # create_frame
else
- create_frame([coords]) # Should not be (coords)
+ if search_mode != 1 || nsearch > 1
+ # The first frame for geometry search has the same coordinates as input
+ create_frame([coords]) # Should not be (coords)
+ end
end
- set_property("energy", energy)
+ set_property("energy", energy) if energy
+ elsif line =~ /BEGINNING GEOMETRY SEARCH POINT/
+ energy = nil # New search has begun, so clear the energy
+ search_mode = 1
+ elsif line =~ /CONSTRAINED OPTIMIZATION POINT/
+ energy = nil # New search has begun, so clear the energy
+ search_mode = 2
elsif line =~ /FINAL .* ENERGY IS *([-.0-9]+) AFTER/
- energy = $1.to_f
- puts energy
+ if search_mode != 2
+ energy = $1.to_f
+ set_property("energy", energy)
+ end
+ elsif line =~ /TOTAL ENERGY += +([-.0-9]+)/
+ energy = $1.to_f
elsif false && line =~ /EQUILIBRIUM GEOMETRY LOCATED/i
set_progress_message(mes + "\nReading optimized coordinates...")
fp.gets; fp.gets; fp.gets
n += 1
break if n >= natoms
end
- if ne_alpha > 0 && ne_beta > 0
- # Allocate basis set record again, to update the atomic coordinates
- allocate_basis_set_record(rflag, ne_alpha, ne_beta)
- end
+ # if ne_alpha > 0 && ne_beta > 0
+ # # Allocate basis set record again, to update the atomic coordinates
+ # allocate_basis_set_record(rflag, ne_alpha, ne_beta)
+ # end
elsif line =~ /ATOMIC BASIS SET/
- while (line = fp.gets)
- break if line =~ /SHELL\s+TYPE\s+PRIMITIVE/
- end
- line = fp.gets
- i = -1
- nprims = 0
- sym = -10 # undefined
- ncomps = 0
+ lines = []
+ lineno = fp.lineno
while (line = fp.gets)
break if line =~ /TOTAL NUMBER OF BASIS SET/
line.chomp!
- if line =~ /^\s*$/
- # End of one shell
- add_gaussian_orbital_shell(sym, nprims, i)
- # puts "add_gaussian_orbital_shell #{sym}, #{nprims}, #{i}"
- nprims = 0
- sym = -10
- next
- end
- a = line.split
- if a.length == 1
- i += 1
- line = fp.gets # Skip the blank line
- next
- elsif a.length == 5 || a.length == 6
- if sym == -10
- case a[1]
- when "S"
- sym = 0; n = 1
- when "P"
- sym = 1; n = 3
- when "L"
- sym = -1; n = 4
- when "D"
- sym = 2; n = 6
- when "F"
- sym = 3; n = 10
- when "G"
- sym = 4; n = 15
- else
- raise MolbyError, "Unknown gaussian shell type at line #{fp.lineno}"
- end
- ncomps += n
- end
- if (a.length == 5 && sym == -1) || (a.length == 6 && sym != -1)
- raise MolbyError, "Wrong format in gaussian shell information at line #{fp.lineno}"
- end
- exp = Float(a[3])
- c = Float(a[4])
- csp = Float(a[5] || 0.0)
- add_gaussian_primitive_coefficients(exp, c, csp)
- nprims += 1
- # puts "add_gaussian_primitive_coefficients #{exp}, #{c}, #{csp}"
- else
- raise MolbyError, "Error in reading basis set information at line #{fp.lineno}"
- end
+ lines.push(line)
end
+ ncomps = sub_load_gamess_log_basis_set(lines, lineno)
elsif line =~ /NUMBER OF OCCUPIED ORBITALS/
line =~ /=\s*(\d+)/
n = Integer($1)
alpha_beta = $1
elsif line =~ /^\s*(EIGENVECTORS|MOLECULAR ORBITALS)\s*$/
if mo_count == 0
- allocate_basis_set_record(rflag, ne_alpha, ne_beta)
- # puts "allocate_basis_set_record #{rflag}, #{ne_alpha}, #{ne_beta}"
+ clear_mo_coefficients
+ set_mo_info(:type=>["UHF", "RHF", "ROHF"][rflag], :alpha=>ne_alpha, :beta=>ne_beta)
end
mo_count += 1
idx = 0
end
n = 0
nf = 0
+ energy = nil
use_input_orientation = false
show_progress_panel("Loading Gaussian out file...")
while 1
if line == nil
break
end
- line.chomp
+ line.chomp!
if line =~ /(Input|Standard) orientation/
match = $1
if match == "Input"
# }
guess_bonds
new_unit = false
- create_frame
+ # create_frame
else
create_frame([coords]) # Should not be (coords)
end
+ if energy
+ # TODO: to ensure whether the energy output line comes before
+ # or after the atomic coordinates.
+ set_property("energy", energy)
+ end
nf += 1
+ elsif line =~ /SCF Done: *E\(\w+\) *= *([-.0-9]+)/
+ energy = $1.to_f
end
end
+ if energy
+ set_property("energy", energy)
+ end
hide_progress_panel
(n > 0 ? true : false)
end
((RubyDialogFrame *)dref)->StartIntervalTimer(-1);
((RubyDialogFrame *)dref)->Show(true);
((RubyDialogFrame *)dref)->Raise();
+#if defined(__WXMAC__)
+ {
+ extern void AddWindowsItemWithTitle(const char *title);
+ wxString str = ((RubyDialogFrame *)dref)->GetLabel();
+ AddWindowsItemWithTitle(str.mb_str(WX_DEFAULT_CONV));
+ }
+#endif
}
void
wxMenuImpl* c = new wxMenuCocoaImpl( peer, menu );
return c;
}
+
+/* 20140614 Toshi Nagata */
+/* Add the frontmost window to the window menu with the given title */
+/* (for RubyDialog) */
+void
+AddWindowsItemWithTitle(const char *title)
+{
+ id win = [NSApp keyWindow];
+ if (win != nil) {
+ [NSApp addWindowsItem:win title:[NSString stringWithUTF8String:title] filename:NO];
+ }
+}
+/* End Toshi Nagata */
OSDN Git Service
