}
/* Allocate ring buffer */
- if (arena->ring != NULL)
+ if (arena->ring != NULL) {
+ free(arena->ring->buf);
free(arena->ring);
- arena->nringframes = 2000 / mol->natoms;
- if (arena->nringframes == 0)
- arena->nringframes = 1;
- arena->ring = (Vector *)malloc(sizeof(Vector) * mol->natoms * arena->nringframes);
+ }
+ arena->ring = (MDRing *)calloc(sizeof(MDRing), 1);
if (arena->ring == NULL)
md_panic(arena, ERROR_out_of_memory);
- arena->ring_next = 0;
- arena->ring_count = 0;
+ arena->ring->size = mol->natoms;
+ if (arena->pressure != NULL && arena->pressure->disabled == 0)
+ arena->ring->size += 4;
+ arena->ring->nframes = 2000 / arena->ring->size;
+ if (arena->ring->nframes < 2)
+ arena->ring->nframes = 2;
+ arena->ring->buf = (Vector *)calloc(sizeof(Vector), arena->ring->size * arena->ring->nframes);
+ if (arena->ring->buf == NULL)
+ md_panic(arena, ERROR_out_of_memory);
+ arena->ring->next = 0;
+ arena->ring->count = 0;
/* Initialize temperature statistics */
arena->sum_temperature = 0.0;
#define get_group_flag(gf, n) (n < (gf)->natoms ? (gf)->flags[(n)/8] & group_flag_mask[(n)%8] : 0)
#define set_group_flag(gf, n, bool) (n < (gf)->natoms ? (bool ? ((gf)->flags[(n)/8] |= group_flag_mask[(n)%8]) : ((gf)->flags[(n)/8] &= ~(group_flag_mask[(n)%8]))) : 0)
+/* Ring buffer for sending coordinates to another thread */
+typedef struct MDRing {
+ Vector *buf; /* Vector buffer */
+ Int size; /* Number of vectors per frame; natoms + (use_cell ? 4 : 0) */
+ Int nframes; /* Number of frames in the ring buffer (2000 / natoms) */
+ Int next; /* Next frame index to store data */
+ Int count; /* Number of frames currently in the ring buffer */
+} MDRing;
+
/* Everything needed for MD */
typedef struct MDArena {
Int refCount;
Int natoms_uniq; /* Number of symmetry-unique atoms */
- Vector *ring; /* Ring buffer for sending coordinates between threads */
- Int nringframes; /* Number of frames in the ring buffer (2000 / natoms) */
- Int ring_next; /* Next frame index to store data */
- Int ring_count; /* Number of frames currently in the ring buffer */
+ MDRing *ring;
/* Parameters are copied from mol->par and gBuiltinParameters for each call to md_prepare() */
Parameter *par;
const char *gMolActionRotateAtoms = "rotate:vdvG";
const char *gMolActionTransformAtoms = "transform:tG";
const char *gMolActionSetAtomPositions = "atomPosition:GV";
-const char *gMolActionInsertFrames = "insertFrames:GV";
+const char *gMolActionInsertFrames = "insertFrames:GVV";
const char *gMolActionRemoveFrames = "removeFrames:G";
const char *gMolActionSetSelection = "selection:G";
const char *gMolActionChangeResidueNumber = "changeResSeq:Gi";
needsSymmetryAmendment = 1;
} else if (strcmp(action->name, gMolActionInsertFrames) == 0) {
int old_nframes, new_nframes;
+ Vector *vp2;
ig = action->args[0].u.igval;
vp = (Vector *)action->args[1].u.arval.ptr;
+ vp2 = (Vector *)action->args[2].u.arval.ptr;
n1 = IntGroupGetCount(ig);
if (n1 == 0)
return 0; /* Do nothing */
if (vp != NULL && action->args[1].u.arval.nitems != n1 * mol->natoms)
return -1; /* Internal inconsistency */
+ if (vp2 != NULL && action->args[2].u.arval.nitems != n1 * 4)
+ return -1; /* Internal inconsistency */
old_nframes = MoleculeGetNumberOfFrames(mol);
- if (MoleculeInsertFrames(mol, ig, vp) < 0)
+ if (MoleculeInsertFrames(mol, ig, vp, vp2) < 0)
return -1; /* Error */
new_nframes = MoleculeGetNumberOfFrames(mol);
if (old_nframes + n1 < new_nframes) {
act2 = MolActionNew(gMolActionRemoveFrames, ig);
act2->frame = mol->cframe;
} else if (strcmp(action->name, gMolActionRemoveFrames) == 0) {
+ Vector *vp2;
ig = action->args[0].u.igval;
n1 = IntGroupGetCount(ig);
if (n1 == 0)
return 0; /* Do nothing */
vp = (Vector *)calloc(sizeof(Vector), n1 * mol->natoms);
- if (MoleculeRemoveFrames(mol, ig, vp) < 0)
+ if (mol->cell != NULL && mol->frame_cells != NULL)
+ vp2 = (Vector *)calloc(sizeof(Vector) * 4, n1);
+ else vp2 = NULL;
+ if (MoleculeRemoveFrames(mol, ig, vp, vp2) < 0)
return -1; /* Error */
- act2 = MolActionNew(gMolActionInsertFrames, ig, n1 * mol->natoms, vp);
+ act2 = MolActionNew(gMolActionInsertFrames, ig, n1 * mol->natoms, vp, (vp2 != NULL ? n1 * 4 : 0), vp2);
act2->frame = mol->cframe;
free(vp);
+ if (vp2 != NULL)
+ free(vp2);
} else if (strcmp(action->name, gMolActionSetSelection) == 0) {
IntGroup *ig2;
ig2 = MoleculeGetSelection(mol);
vbuf[i].z = dval[3];
}
ig = IntGroupNewWithPoints(nframes, 1, -1);
- MolActionCreateAndPerform(mol, gMolActionInsertFrames, ig, natoms, vbuf);
+ MolActionCreateAndPerform(mol, gMolActionInsertFrames, ig, natoms, vbuf, 0, NULL);
IntGroupRelease(ig);
nframes++;
if (n1) {
vpp->z = zp[i];
}
}
- if (MoleculeInsertFrames(mp, ig, vp) < 0)
+ /* TODO: implement frame-specific cells */
+ if (MoleculeInsertFrames(mp, ig, vp, NULL) < 0)
snprintf(errbuf, errbufsize, "Cannot insert frames");
if (dcd.with_unitcell) {
Vector ax, ay, az, orig;
return -1; /* Non-regular transform */
/* Calculate the reciprocal cell parameters */
- cp->rcell[0] = sqrt(cp->rtr[0] * cp->rtr[0] + cp->rtr[1] * cp->rtr[1] + cp->rtr[2] * cp->rtr[2]);
- cp->rcell[1] = sqrt(cp->rtr[3] * cp->rtr[3] + cp->rtr[4] * cp->rtr[4] + cp->rtr[5] * cp->rtr[5]);
- cp->rcell[2] = sqrt(cp->rtr[6] * cp->rtr[6] + cp->rtr[7] * cp->rtr[7] + cp->rtr[8] * cp->rtr[8]);
- cp->rcell[3] = acos((cp->rtr[3] * cp->rtr[6] + cp->rtr[4] * cp->rtr[7] + cp->rtr[5] * cp->rtr[8]) / (cp->rcell[1] * cp->rcell[2])) * kRad2Deg;
- cp->rcell[4] = acos((cp->rtr[6] * cp->rtr[0] + cp->rtr[7] * cp->rtr[1] + cp->rtr[8] * cp->rtr[2]) / (cp->rcell[2] * cp->rcell[0])) * kRad2Deg;
- cp->rcell[5] = acos((cp->rtr[0] * cp->rtr[3] + cp->rtr[1] * cp->rtr[4] + cp->rtr[2] * cp->rtr[5]) / (cp->rcell[0] * cp->rcell[1])) * kRad2Deg;
+ cp->rcell[0] = sqrt(cp->rtr[0] * cp->rtr[0] + cp->rtr[3] * cp->rtr[3] + cp->rtr[6] * cp->rtr[6]);
+ cp->rcell[1] = sqrt(cp->rtr[1] * cp->rtr[1] + cp->rtr[4] * cp->rtr[4] + cp->rtr[7] * cp->rtr[7]);
+ cp->rcell[2] = sqrt(cp->rtr[2] * cp->rtr[2] + cp->rtr[5] * cp->rtr[5] + cp->rtr[8] * cp->rtr[8]);
+ cp->rcell[3] = acos((cp->rtr[1] * cp->rtr[2] + cp->rtr[4] * cp->rtr[5] + cp->rtr[7] * cp->rtr[8]) / (cp->rcell[1] * cp->rcell[2])) * kRad2Deg;
+ cp->rcell[4] = acos((cp->rtr[2] * cp->rtr[0] + cp->rtr[5] * cp->rtr[3] + cp->rtr[8] * cp->rtr[6]) / (cp->rcell[2] * cp->rcell[0])) * kRad2Deg;
+ cp->rcell[5] = acos((cp->rtr[0] * cp->rtr[1] + cp->rtr[3] * cp->rtr[4] + cp->rtr[6] * cp->rtr[7]) / (cp->rcell[0] * cp->rcell[1])) * kRad2Deg;
- if (!calc_abc) {
+ if (calc_abc) {
/* Calculate a, b, c, alpha, beta, gamma */
- cp->cell[0] = sqrt(cp->tr[0] * cp->tr[0] + cp->tr[1] * cp->tr[1] + cp->tr[2] * cp->tr[2]);
- cp->cell[1] = sqrt(cp->tr[3] * cp->tr[3] + cp->tr[4] * cp->tr[4] + cp->tr[5] * cp->tr[5]);
- cp->cell[2] = sqrt(cp->tr[6] * cp->tr[6] + cp->tr[7] * cp->tr[7] + cp->tr[8] * cp->tr[8]);
- cp->cell[3] = acos((cp->tr[3] * cp->tr[6] + cp->tr[4] * cp->tr[7] + cp->tr[5] * cp->tr[8]) / (cp->cell[1] * cp->cell[2])) * kRad2Deg;
- cp->cell[4] = acos((cp->tr[6] * cp->tr[0] + cp->tr[7] * cp->tr[1] + cp->tr[8] * cp->tr[2]) / (cp->cell[2] * cp->cell[0])) * kRad2Deg;
- cp->cell[5] = acos((cp->tr[0] * cp->tr[3] + cp->tr[1] * cp->tr[4] + cp->tr[2] * cp->tr[5]) / (cp->cell[0] * cp->cell[1])) * kRad2Deg;
+ cp->cell[0] = sqrt(cp->tr[0] * cp->tr[0] + cp->tr[3] * cp->tr[3] + cp->tr[6] * cp->tr[6]);
+ cp->cell[1] = sqrt(cp->tr[1] * cp->tr[1] + cp->tr[4] * cp->tr[4] + cp->tr[7] * cp->tr[7]);
+ cp->cell[2] = sqrt(cp->tr[2] * cp->tr[2] + cp->tr[5] * cp->tr[5] + cp->tr[8] * cp->tr[8]);
+ cp->cell[3] = acos((cp->tr[1] * cp->tr[2] + cp->tr[4] * cp->tr[5] + cp->tr[7] * cp->tr[8]) / (cp->cell[1] * cp->cell[2])) * kRad2Deg;
+ cp->cell[4] = acos((cp->tr[2] * cp->tr[0] + cp->tr[5] * cp->tr[3] + cp->tr[8] * cp->tr[6]) / (cp->cell[2] * cp->cell[0])) * kRad2Deg;
+ cp->cell[5] = acos((cp->tr[0] * cp->tr[1] + cp->tr[3] * cp->tr[4] + cp->tr[6] * cp->tr[7]) / (cp->cell[0] * cp->cell[1])) * kRad2Deg;
}
return 0;
}
int
-MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame)
+MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame, const Vector *inFrameCell)
{
- int i, j, count, n_new, n_old, natoms, exframes, last_inserted;
+ int i, j, count, n_new, n_old, natoms, exframes, last_inserted, old_count;
Vector *tempv, *vp;
Atom *ap;
if (mp == NULL || (natoms = mp->natoms) == 0 || (count = IntGroupGetCount(group)) <= 0)
n_new += exframes;
} else exframes = 0;
- tempv = (Vector *)malloc(sizeof(Vector) * n_new);
+ tempv = (Vector *)malloc(sizeof(Vector) * n_new * 4); /* "*4" for handling cells */
if (tempv == NULL)
return -1;
vp[j] = ap->r;
ap->frames = vp;
}
+ if (mp->cell != NULL && (mp->frame_cells != NULL || inFrameCell != NULL)) {
+ if (mp->frame_cells == NULL)
+ vp = (Vector *)calloc(sizeof(Vector), n_new * 4);
+ else
+ vp = (Vector *)realloc(mp->frame_cells, sizeof(Vector) * 4 * n_new);
+ if (vp == NULL) {
+ __MoleculeUnlock(mp);
+ return -1;
+ }
+ mp->frame_cells = vp;
+ }
/* group = [n0..n1-1, n2..n3-1, ...] */
/* s = t = 0, */
...
tempv[nl..n_new-1] <- ap[s..s+(n_new-nl-1)], s += n_new-nl
At last, s will become n_old and t will become count. */
- for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
- int s, t, ns, ne;
- Vector cr = ap->r;
+ for (i = 0, ap = mp->atoms; i <= mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ int s, t, ns, ne, mult;
+ Vector cr;
ne = s = t = 0;
- vp = ap->frames;
+ if (i == mp->natoms) {
+ if (mp->cell == NULL || mp->frame_cells == NULL)
+ break;
+ vp = mp->frame_cells;
+ mult = 4;
+ } else {
+ cr = ap->r;
+ vp = ap->frames;
+ mult = 1;
+ }
for (j = 0; (ns = IntGroupGetStartPoint(group, j)) >= 0; j++) {
if (ns > ne) {
- memmove(tempv + ne, vp + s, sizeof(Vector) * (ns - ne));
+ memmove(tempv + ne * mult, vp + s * mult, sizeof(Vector) * mult * (ns - ne));
s += ns - ne;
}
ne = IntGroupGetEndPoint(group, j);
while (ns < ne) {
- if (inFrame != NULL)
- tempv[ns] = inFrame[natoms * t + i];
- else tempv[ns] = cr;
+ if (i == mp->natoms) {
+ if (inFrameCell != NULL) {
+ tempv[ns * 4] = inFrameCell[t * 4];
+ tempv[ns * 4 + 1] = inFrameCell[t * 4 + 1];
+ tempv[ns * 4 + 2] = inFrameCell[t * 4 + 2];
+ tempv[ns * 4 + 3] = inFrameCell[t * 4 + 3];
+ } else {
+ tempv[ns * 4] = mp->cell->axes[0];
+ tempv[ns * 4 + 1] = mp->cell->axes[1];
+ tempv[ns * 4 + 2] = mp->cell->axes[2];
+ tempv[ns * 4 + 3] = mp->cell->origin;
+ }
+ } else {
+ if (inFrame != NULL)
+ tempv[ns] = inFrame[natoms * t + i];
+ else
+ tempv[ns] = cr;
+ }
t++;
ns++;
}
}
if (n_new > ne) {
- memmove(tempv + ne, vp + s, sizeof(Vector) * (n_new - ne));
+ memmove(tempv + ne * mult, vp + s * mult, sizeof(Vector) * mult * (n_new - ne));
s += n_new - ne;
}
- ap->nframes = n_new;
- memmove(vp, tempv, sizeof(Vector) * n_new);
+ if (i < mp->natoms)
+ ap->nframes = n_new;
+ memmove(vp, tempv, sizeof(Vector) * mult * n_new);
}
free(tempv);
mp->nframes = n_new;
}
int
-MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame)
+MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame, Vector *outFrameCell)
{
- int i, count, n_new, n_old, natoms, nframes;
+ int i, count, n_new, n_old, natoms, nframes, old_count;
Vector *tempv, *vp;
Atom *ap;
if (mp == NULL || (natoms = mp->natoms) == 0 || (count = IntGroupGetCount(group)) <= 0)
/* outFrame[] should have enough size for Vector * natoms * group.count */
memset(outFrame, 0, sizeof(Vector) * natoms * count);
+ if (mp->cell != NULL && mp->frame_cells != NULL)
+ memset(outFrameCell, 0, sizeof(Vector) * 4 * count);
n_old = MoleculeGetNumberOfFrames(mp);
n_new = n_old - count;
if (n_new < 0)
n_new = 0;
- tempv = (Vector *)malloc(sizeof(Vector) * n_old);
+ tempv = (Vector *)malloc(sizeof(Vector) * n_old * 4); /* "*4" for handling cells */
if (tempv == NULL)
return -1;
At last, s will become n_new and t will become count. */
__MoleculeLock(mp);
nframes = 0;
- for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ for (i = 0, ap = mp->atoms; i <= mp->natoms; i++, ap = ATOM_NEXT(ap)) {
int s, t, j, ns, ne;
- /* Copy ap->frames to tempv */
- memset(tempv, 0, sizeof(Vector) * n_old);
- vp = ap->frames;
- if (vp != NULL)
- memmove(tempv, vp, sizeof(Vector) * (ap->nframes > n_old ? n_old : ap->nframes));
- else {
- ap->frames = vp = (Vector *)calloc(sizeof(Vector), n_old);
+ int mult;
+ /* if i == mp->natoms, mp->frame_cells is handled */
+ if (i == mp->natoms) {
+ if (mp->cell == NULL || mp->frame_cells == NULL)
+ break;
+ mult = 4;
+ vp = mp->frame_cells;
+ old_count = n_old;
+ } else {
+ mult = 1;
+ vp = ap->frames;
if (vp == NULL) {
- __MoleculeUnlock(mp);
- return -1;
+ ap->frames = vp = (Vector *)calloc(sizeof(Vector), n_old);
+ if (vp == NULL) {
+ __MoleculeUnlock(mp);
+ return -1;
+ }
}
+ old_count = ap->nframes;
}
+
+ /* Copy vp to tempv */
+ memset(tempv, 0, sizeof(Vector) * mult * n_old);
+ memmove(tempv, vp, sizeof(Vector) * mult * (old_count > n_old ? n_old : old_count));
ne = ns = s = t = 0;
for (j = 0; ns < n_old && (ns = IntGroupGetStartPoint(group, j)) >= 0; j++) {
if (ns > n_old)
ns = n_old;
if (ns > ne) {
- memmove(vp + s, tempv + ne, sizeof(Vector) * (ns - ne));
+ memmove(vp + s * mult, tempv + ne * mult, sizeof(Vector) * mult * (ns - ne));
s += ns - ne;
}
ne = IntGroupGetEndPoint(group, j);
if (ne > n_old)
ne = n_old;
while (ns < ne) {
- outFrame[natoms * t + i] = tempv[ns];
+ if (i < mp->natoms)
+ outFrame[natoms * t + i] = tempv[ns];
+ else if (outFrameCell != NULL) {
+ outFrameCell[i * 4] = tempv[ns * 4];
+ outFrameCell[i * 4 + 1] = tempv[ns * 4 + 1];
+ outFrameCell[i * 4 + 2] = tempv[ns * 4 + 2];
+ outFrameCell[i * 4 + 3] = tempv[ns * 4 + 3];
+ }
t++;
ns++;
}
}
if (n_old > ne) {
- memmove(vp + s, tempv + ne, sizeof(Vector) * (n_old - ne));
+ memmove(vp + s * mult, tempv + ne * mult, sizeof(Vector) * mult * (n_old - ne));
s += n_old - ne;
}
- ap->nframes = s;
+ if (i < mp->natoms)
+ ap->nframes = s;
if (nframes < s)
nframes = s;
if (s == 0) {
- free(ap->frames);
- ap->frames = NULL;
+ if (i < mp->natoms) {
+ free(ap->frames);
+ ap->frames = NULL;
+ } else {
+ free(mp->frame_cells);
+ mp->frame_cells = NULL;
+ }
} else {
- ap->frames = (Vector *)realloc(ap->frames, sizeof(Vector) * s);
+ if (i < mp->natoms)
+ ap->frames = (Vector *)realloc(ap->frames, sizeof(Vector) * s);
+ else
+ mp->frame_cells = (Vector *)realloc(mp->frame_cells, sizeof(Vector) * 4 * s);
}
}
free(tempv);
return count;
}
+#if 0
int
MoleculeInsertFrame(Molecule *mp, int index, const Vector *inFrame)
{
__MoleculeUnlock(mp);
return frame;
}
+#endif
int
MoleculeSelectFrame(Molecule *mp, int frame, int copyback)
ok = 1;
}
}
+
+ if (mp->cell != NULL && mp->frame_cells != NULL) {
+ 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);
+ }
+
__MoleculeUnlock(mp);
if (ok) {
mp->cframe = frame;
Int nframes; /* The number of frames. This is a cached value, and should be
recalculated from the atoms if it is -1 */
Int cframe; /* The current frame number */
+
+ Vector *frame_cells; /* The cell vectors for frames; (nframes*3) array of Vectors */
+
struct MainView *mview; /* Reference to the MainView object if present (no retain) */
Int modifyCount; /* Internal counter for modification. This value is not to be modified
manually; instead, call MoleculeIncrementModifyCount() whenever
int MoleculeIsFragmentRotatable(Molecule *mp, IntGroup *group, int *n1, int *n2, IntGroup **rotGroup);
int MoleculeGetNumberOfFrames(Molecule *mp);
-int MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame);
-int MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame);
-int MoleculeInsertFrame(Molecule *mp, int index, const Vector *inFrame);
-int MoleculeRemoveFrame(Molecule *mp, int frame, Vector *outFrame);
+int MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame, const Vector *inFrameCell);
+int MoleculeRemoveFrames(Molecule *mp, IntGroup *group, Vector *outFrame, Vector *outFrameCell);
+/*int MoleculeInsertFrame(Molecule *mp, int index, const Vector *inFrame);
+int MoleculeRemoveFrame(Molecule *mp, int frame, Vector *outFrame); */
int MoleculeSelectFrame(Molecule *mp, int frame, int copyback);
int MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, const Vector *dyp, const Vector *dzp, Int nx, Int ny, Int nz, int (*callback)(double progress, void *ref), void *ref);
/*
* call-seq:
- * insert_frame(integer, coordinates = nil) -> bool
- * insert_frames(intGroup = nil, coordinates = nil) -> bool
+ * 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
static VALUE
s_Molecule_InsertFrames(int argc, VALUE *argv, VALUE self)
{
- VALUE val, coords;
+ VALUE val, coords, cells;
Molecule *mol;
IntGroup *ig;
int count, ival, i, j, len, nframes;
- Vector *vp;
+ Vector *vp, *vp2;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "11", &val, &coords);
+ 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");
+ if (len != RARRAY_LEN(coords))
+ rb_raise(rb_eTypeError, "the number of entries for coordinates and cell axes does not match");
+ }
nframes = MoleculeGetNumberOfFrames(mol);
if (val == Qnil) {
ig = IntGroupNewWithPoints(nframes, (len > 0 ? len : 1), -1);
}
count = IntGroupGetCount(ig);
vp = ALLOC_N(Vector, mol->natoms * count);
+ if (cells != Qnil && len > 0)
+ vp2 = ALLOC_N(Vector, 4 * count);
+ else vp2 = NULL;
if (len > 0) {
- VALUE *ptr;
+ VALUE *ptr, *ptr2;
+ int j, len2;
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 < len; i++) {
- VALUE *ptr2;
- int j, len2;
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]);
for (j = 0; j < mol->natoms; j++)
VectorFromValue(ptr2[j], &vp[i * mol->natoms + j]);
}
+ if (vp2 != NULL) {
+ ptr = RARRAY_PTR(cells);
+ for (i = 0; i < len; 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]);
+ }
+ }
} else {
Atom *ap;
for (i = 0; i < count; i++) {
}
}
}
- ival = MolActionCreateAndPerform(mol, gMolActionInsertFrames, ig, mol->natoms * count, vp);
+ ival = MolActionCreateAndPerform(mol, gMolActionInsertFrames, ig, mol->natoms * count, vp, (vp2 != NULL ? 4 * count : 0), vp2);
IntGroupRelease(ig);
free(vp);
+ if (vp2 != NULL)
+ free(vp2);
return (ival >= 0 ? val : Qnil);
}
if (arena->step > start_step) {
/* Create a new frame and copy new coordinates */
ig = IntGroupNewWithPoints(MoleculeGetNumberOfFrames(arena->xmol), 1, -1);
- MolActionCreateAndPerform(arena->xmol, gMolActionInsertFrames, ig, 0, NULL);
+ MolActionCreateAndPerform(arena->xmol, gMolActionInsertFrames, ig, 0, NULL, 0, NULL);
IntGroupRelease(ig);
md_copy_coordinates_from_internal(arena);
}
DialogItem は C レベルのデータは持たず、単なる cDialogItem のインスタンスとする。属性は、今まで Dialog 中のハッシュに保持していたが、DialogItem のインスタンス変数として実装する。
Dialog#item(:type, hash) でアイテムを作成。Dialog#set_attr, Dialog#get_attr はやめて、DialogItem#[], DialogItem#[]= で実装する。DialogItem#属性名, DialogItem#属性名= も実装する（ParameterRef と同様）。また、Dialog#action の引数はアイテム番号でなく DialogItem オブジェクトとする。
DialogItem の action 属性は今までと同様だが、引数としてはアイテム番号でなく DialogItem オブジェクトそのものが渡される。
-
\ No newline at end of file
+
+2010.3.21.
+ Pressure control を使えるようにしたい。フレームごとに unit cell パラメータを保持できるようにしないといけない。
if (mol->requestAbortThread)
r = -1;
else {
-
/* Copy the coordinate to the ring buffer */
- Vector *rp = mol->arena->ring + mol->natoms * mol->arena->ring_next;
+ MDRing *ring = mol->arena->ring;
+ Vector *rp = ring->buf + ring->size * ring->next;
Int j;
Atom *ap;
MoleculeLock(mol);
for (j = 0, ap = mol->arena->mol->atoms; j < mol->natoms; j++, ap = ATOM_NEXT(ap)) {
rp[j] = ap->r;
}
- mol->arena->ring_next = (mol->arena->ring_next + 1) % mol->arena->nringframes;
- if (mol->arena->ring_count < mol->arena->nringframes)
- mol->arena->ring_count++;
+ if (j < ring->size) {
+ XtalCell *cp = mol->arena->mol->cell;
+ if (cp != NULL) {
+ rp[j++] = cp->axes[0];
+ rp[j++] = cp->axes[1];
+ rp[j++] = cp->axes[2];
+ rp[j++] = cp->origin;
+ }
+ }
+ ring->next = (ring->next + 1) % ring->nframes;
+ if (ring->count < ring->nframes)
+ ring->count++;
MoleculeUnlock(mol);
/* Create a new frame and copy the new coordinates */
{
Int i, j, n, old_nframes;
Atom *ap;
+ MDRing *ring;
/* Create new frame(s) and copy the new coordinates from the ring buffer */
MoleculeLock(mol);
- n = mol->arena->ring_count;
+ ring = mol->arena->ring;
+ n = ring->count;
if (n > 0) {
IntGroup *ig;
Vector *rp;
old_nframes = MoleculeGetNumberOfFrames(mol);
+ /* It is more convenient to set cell parameter when inserting frames, whereas
+ the coordinates can be set afterwards */
+ if (ring->size > mol->natoms) {
+ rp = (Vector *)calloc(sizeof(Vector) * 4, n);
+ for (i = 0; i < n; i++) {
+ j = ((ring->next - n + i + ring->nframes) % ring->nframes) * ring->size + mol->natoms;
+ rp[i * 4] = ring->buf[j++];
+ rp[i * 4 + 1] = ring->buf[j++];
+ rp[i * 4 + 2] = ring->buf[j++];
+ rp[i * 4 + 3] = ring->buf[j++];
+ }
+ } else rp = NULL;
ig = IntGroupNewWithPoints(old_nframes, n, -1);
- MolActionCreateAndPerform(mol, gMolActionInsertFrames, ig, 0, NULL);
+ MolActionCreateAndPerform(mol, gMolActionInsertFrames, ig, 0, NULL, (rp != NULL ? n * 4 : 0), rp);
+ if (rp != NULL)
+ free(rp);
IntGroupRelease(ig);
for (i = 0; i < n; i++) {
MoleculeSelectFrame(mol, old_nframes + i, 1);
- rp = mol->arena->ring + ((mol->arena->ring_next + mol->arena->nringframes - n + i) % mol->arena->nringframes) * mol->natoms;
+ rp = ring->buf + ((ring->next - n + i + ring->nframes) % ring->nframes) * ring->size;
for (j = 0, ap = mol->atoms; j < mol->natoms; j++, ap = ATOM_NEXT(ap))
ap->r = rp[j];
}
- mol->arena->ring_count = 0;
+ ring->count = 0;
}
MoleculeUnlock(mol);
}
OSDN Git Service
