mol->syms = NULL;
}
} else if (strcmp(action->name, gMolActionSetCell) == 0) {
- double *dp, d[6];
- if (mol->cell == NULL)
+ double *dp, d[12];
+ int convertCoord, n2;
+ if (mol->cell == NULL) {
d[0] = 0.0;
- else {
+ n1 = 0;
+ } else {
for (n1 = 0; n1 < 6; n1++)
d[n1] = mol->cell->cell[n1];
+ if (mol->cell->has_sigma) {
+ for (n1 = 6; n1 < 12; n1++)
+ d[n1] = mol->cell->cellsigma[n1 - 6];
+ }
}
- n1 = action->args[1].u.ival;
+ convertCoord = action->args[1].u.ival;
dp = action->args[0].u.arval.ptr;
- if (action->args[0].u.arval.nitems == 0)
- MoleculeSetCell(mol, 0, 0, 0, 0, 0, 0, n1);
- else
- MoleculeSetCell(mol, dp[0], dp[1], dp[2], dp[3], dp[4], dp[5], n1);
- act2 = MolActionNew(gMolActionSetCell, (d[0] == 0.0 ? 0 : 6), d, n1);
+ n2 = action->args[0].u.arval.nitems;
+ if (n2 == 0)
+ MoleculeSetCell(mol, 0, 0, 0, 0, 0, 0, convertCoord);
+ else {
+ MoleculeSetCell(mol, dp[0], dp[1], dp[2], dp[3], dp[4], dp[5], convertCoord);
+ if (n2 == 12) {
+ mol->cell->has_sigma = 1;
+ for (n2 = 6; n2 < 12; n2++)
+ mol->cell->cellsigma[n2 - 6] = dp[n2];
+ } else mol->cell->has_sigma = 0;
+ }
+ act2 = MolActionNew(gMolActionSetCell, n1, (n1 == 0 ? NULL : d), convertCoord);
needsRebuildMDArena = 1;
} else if (strcmp(action->name, gMolActionSetBox) == 0) {
Vector v[4];
atomType = fbuf[6];
if ((atomType >= 0 && atomType <= 5) || (atomType >= 8 && atomType <= 10)) {
/* Anisotropic thermal parameters */
- MoleculeSetAniso(mp, atomIndex, atomType, fbuf[0], fbuf[1], fbuf[2], fbuf[3], fbuf[5], fbuf[4]);
+ MoleculeSetAniso(mp, atomIndex, atomType, fbuf[0], fbuf[1], fbuf[2], fbuf[3], fbuf[5], fbuf[4], NULL);
}
if (ibuf[0] != 0)
section = 3;
Transform tr;
int ibuf[12];
float fbuf[12];
+ Double dbuf[12];
Int nsfacs = 0;
Int nbonds, *bonds;
char (*sfacs)[4] = NULL;
ElementToString(ap->atomicNumber, ap->element); */
}
guessElement(ap);
+ if (n == 10 || fbuf[4] >= 0.0) {
+ int i, c, j;
+ /* Read in the standard deviations */
+ ReadLine(buf, sizeof buf, fp, &lineNumber);
+ for (i = 0; i < 9; i++) {
+ j = 3 + i * 8;
+ c = buf[j + 8];
+ buf[j + 8] = 0;
+ dbuf[i] = strtod(buf + j, NULL);
+ buf[j + 8] = c;
+ }
+ ap->sigma.x = dbuf[0];
+ ap->sigma.y = dbuf[1];
+ ap->sigma.z = dbuf[2];
+ }
if (n == 5)
ap->tempFactor = fbuf[4] * 78.9568352087147; /* 8*pi*pi */
else
- MoleculeSetAniso(mp, atomIndex, 8, fbuf[4], fbuf[5], fbuf[6], fbuf[9], fbuf[7], fbuf[8]);
+ MoleculeSetAniso(mp, atomIndex, 8, fbuf[4], fbuf[5], fbuf[6], fbuf[9], fbuf[7], fbuf[8], dbuf);
ap->resSeq = currentResSeq;
strncpy(ap->resName, currentResName, 4);
}
} else {
cp = mp->cell;
if (cp == NULL) {
- cp = (XtalCell *)malloc(sizeof(XtalCell));
+ cp = (XtalCell *)calloc(sizeof(XtalCell), 1);
if (cp == NULL)
Panic("Low memory during setting cell parameters");
mp->cell = cp;
for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
Aniso *anp = ap->aniso;
if (anp != NULL) {
- MoleculeSetAniso(mp, i, 0, anp->bij[0], anp->bij[1], anp->bij[2], anp->bij[3], anp->bij[4], anp->bij[5]);
+ MoleculeSetAniso(mp, i, 0, anp->bij[0], anp->bij[1], anp->bij[2], anp->bij[3], anp->bij[4], anp->bij[5], anp->bsig);
}
}
__MoleculeUnlock(mp);
}
void
-MoleculeSetAniso(Molecule *mp, int n1, int type, Double x11, Double x22, Double x33, Double x12, Double x23, Double x31)
+MoleculeSetAniso(Molecule *mp, int n1, int type, Double x11, Double x22, Double x33, Double x12, Double x23, Double x31, const Double *sigmaptr)
{
Double d, dx;
int u = 0;
anp = mp->atoms[n1].aniso;
__MoleculeLock(mp);
if (anp == NULL) {
- anp = (Aniso *)malloc(sizeof(Aniso));
+ anp = (Aniso *)calloc(sizeof(Aniso), 1);
if (anp == NULL) {
__MoleculeUnlock(mp);
Panic("Low memory during setting anisotropic atom parameters");
anp->bij[3] = x12 * dx;
anp->bij[4] = x23 * dx;
anp->bij[5] = x31 * dx;
+ if (sigmaptr != NULL) {
+ anp->has_bsig = 1;
+ anp->bsig[0] = sigmaptr[0] * d;
+ anp->bsig[1] = sigmaptr[1] * d;
+ anp->bsig[2] = sigmaptr[2] * d;
+ anp->bsig[3] = sigmaptr[3] * dx;
+ anp->bsig[4] = sigmaptr[4] * dx;
+ anp->bsig[5] = sigmaptr[5] * dx;
+ } else {
+ anp->has_bsig = 0;
+ anp->bsig[0] = anp->bsig[1] = anp->bsig[2] = anp->bsig[3] = anp->bsig[4] = anp->bsig[5] = 0.0;
+ }
cp = mp->cell;
if (cp != NULL && u == 1) {
anp->bij[0] *= cp->rcell[0] * cp->rcell[0];
anp->bij[3] *= cp->rcell[0] * cp->rcell[1]; /* * cos(cp->rcell[5] * kDeg2Rad); */
anp->bij[4] *= cp->rcell[1] * cp->rcell[2]; /* * cos(cp->rcell[3] * kDeg2Rad); */
anp->bij[5] *= cp->rcell[2] * cp->rcell[0]; /* * cos(cp->rcell[4] * kDeg2Rad); */
+ if (sigmaptr != NULL) {
+ anp->bsig[0] *= cp->rcell[0] * cp->rcell[0];
+ anp->bsig[1] *= cp->rcell[1] * cp->rcell[1];
+ anp->bsig[2] *= cp->rcell[2] * cp->rcell[2];
+ anp->bsig[3] *= cp->rcell[0] * cp->rcell[1];
+ anp->bsig[4] *= cp->rcell[1] * cp->rcell[2];
+ anp->bsig[5] *= cp->rcell[2] * cp->rcell[0];
+ }
}
/* Calculate the principal axes (in Cartesian coordinates) */
/* mp->is_xtal_coord = 0; */
return 0;
}
+ memset(&b, 0, sizeof(b));
b.axes[0] = (ax != NULL ? *ax : zeroVec);
b.axes[1] = (ay != NULL ? *ay : zeroVec);
b.axes[2] = (az != NULL ? *az : zeroVec);
__MoleculeLock(mp);
if (mp->cell != NULL)
free(mp->cell);
- mp->cell = (XtalCell *)malloc(sizeof(XtalCell));
+ mp->cell = (XtalCell *)calloc(sizeof(XtalCell), 1);
if (mp->cell != NULL) {
memmove(mp->cell, &b, sizeof(XtalCell));
n = 0;
/* Anisotropic thermal parameter */
typedef struct Aniso {
Double bij[6]; /* b11, b22, b33, b12, b23, b31 (ORTEP type 0) */
+ char has_bsig; /* Has sigma values? */
+ Double bsig[6]; /* sigma values */
Mat33 pmat; /* A 3x3 matrix whose three row vectors are the principal axes of the ellipsoid. Note: If the B matrix is not positive definite, the axis length corresponding to the negative eigenvalue is replaced with 0.001. */
} Aniso;
Vector r; /* position */
Vector v; /* velocity */
Vector f; /* force */
+ Vector sigma; /* For crystallographic data only; sigma for each crystallographic coordinates */
+ /* (Unlike r, these are not converted to the cartesian system) */
Double occupancy;
Double tempFactor;
Aniso *aniso;
Vector axes[3]; /* Cartesian unit vectors along the three axis */
Vector origin; /* Cartesian origin of the periodic box */
char flags[3]; /* 1 for periodic, 0 for non-periodic */
+ char has_sigma; /* Has sigma? */
Transform tr; /* Crystal coord -> cartesian */
Transform rtr; /* Cartesian -> crystal coord */
+ Double cellsigma[6]; /* For crystallographic data; sigma for the cell parameters */
} XtalCell;
/* Periodic box parameter */
void AtomRefRelease(AtomRef *aref);
void MoleculeSetCell(Molecule *mp, Double a, Double b, Double c, Double alpha, Double beta, Double gamma, int convertCoordinates);
-void MoleculeSetAniso(Molecule *mp, int n1, int type, Double x11, Double x22, Double x33, Double x12, Double x23, Double x31);
+void MoleculeSetAniso(Molecule *mp, int n1, int type, Double x11, Double x22, Double x33, Double x12, Double x23, Double x31, const Double *sigmap);
int MoleculeSetPeriodicBox(Molecule *mp, const Vector *ax, const Vector *ay, const Vector *az, const Vector *ao, const char *periodic);
int MoleculeReadCoordinatesFromFile(Molecule *mp, const char *fname, const char *ftype, char *errbuf, int errbufsize);
s_WeightSym, s_ElementSym, s_AtomicNumberSym, s_ConnectsSym,
s_RSym, s_XSym, s_YSym, s_ZSym,
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;
return rb_float_new(s_AtomRef_GetFractionalRAsVector(self).z);
}
+static VALUE s_AtomRef_GetSigma(VALUE self) {
+ return ValueFromVector(&(s_AtomFromValue(self)->sigma));
+}
+
+static VALUE s_AtomRef_GetSigmaX(VALUE self) {
+ return rb_float_new(s_AtomFromValue(self)->sigma.x);
+}
+
+static VALUE s_AtomRef_GetSigmaY(VALUE self) {
+ return rb_float_new(s_AtomFromValue(self)->sigma.y);
+}
+
+static VALUE s_AtomRef_GetSigmaZ(VALUE self) {
+ return rb_float_new(s_AtomFromValue(self)->sigma.z);
+}
+
static VALUE s_AtomRef_GetV(VALUE self) {
return ValueFromVector(&(s_AtomFromValue(self)->v));
}
retval = rb_ary_new();
for (i = 0; i < 6; i++)
rb_ary_push(retval, rb_float_new(ap->aniso->bij[i]));
+ if (ap->aniso->has_bsig) {
+ rb_ary_push(retval, INT2NUM(0));
+ for (i = 0; i < 6; i++)
+ rb_ary_push(retval, rb_float_new(ap->aniso->bsig[i]));
+ }
return retval;
}
return val;
}
+static VALUE s_AtomRef_SetSigma(VALUE self, VALUE val) {
+ Vector v;
+ Molecule *mp;
+ VALUE oval = s_AtomRef_GetSigma(self);
+ VectorFromValue(val, &v);
+ s_AtomAndMoleculeFromValue(self, &mp)->sigma = v;
+ s_RegisterUndoForAtomAttrChange(self, s_SigmaSym, val, oval);
+ mp->needsMDCopyCoordinates = 1;
+ return val;
+}
+
+static VALUE s_AtomRef_SetSigmaX(VALUE self, VALUE val) {
+ Double f;
+ Molecule *mp;
+ VALUE oval = s_AtomRef_GetSigmaX(self);
+ val = rb_Float(val);
+ f = NUM2DBL(val);
+ s_AtomAndMoleculeFromValue(self, &mp)->sigma.x = f;
+ s_RegisterUndoForAtomAttrChange(self, s_SigmaXSym, val, oval);
+ mp->needsMDCopyCoordinates = 1;
+ return val;
+}
+
+static VALUE s_AtomRef_SetSigmaY(VALUE self, VALUE val) {
+ Double f;
+ Molecule *mp;
+ VALUE oval = s_AtomRef_GetSigmaY(self);
+ val = rb_Float(val);
+ f = NUM2DBL(val);
+ s_AtomAndMoleculeFromValue(self, &mp)->sigma.y = f;
+ s_RegisterUndoForAtomAttrChange(self, s_SigmaYSym, val, oval);
+ mp->needsMDCopyCoordinates = 1;
+ return val;
+}
+
+static VALUE s_AtomRef_SetSigmaZ(VALUE self, VALUE val) {
+ Double f;
+ Molecule *mp;
+ VALUE oval = s_AtomRef_GetSigmaZ(self);
+ val = rb_Float(val);
+ f = NUM2DBL(val);
+ s_AtomAndMoleculeFromValue(self, &mp)->sigma.z = f;
+ s_RegisterUndoForAtomAttrChange(self, s_SigmaZSym, val, oval);
+ mp->needsMDCopyCoordinates = 1;
+ return val;
+}
static VALUE s_AtomRef_SetV(VALUE self, VALUE val) {
Vector v;
AtomRef *aref;
int i, n, type;
VALUE *valp;
- Double f[6];
+ Double f[12];
VALUE oval = s_AtomRef_GetAniso(self);
Data_Get_Struct(self, AtomRef, aref);
val = rb_funcall(val, rb_intern("to_a"), 0);
if (n >= 7)
type = NUM2INT(rb_Integer(valp[6]));
else type = 0;
+ if (n >= 13) {
+ for (i = 0; i < 6; i++)
+ f[i + 6] = NUM2DBL(rb_Float(valp[i + 7]));
+ } else {
+ for (i = 0; i < 6; i++)
+ f[i + 6] = 0.0;
+ }
i = s_AtomIndexFromValue(self, NULL, NULL);
- MoleculeSetAniso(aref->mol, i, type, f[0], f[1], f[2], f[3], f[4], f[5]);
+ MoleculeSetAniso(aref->mol, i, type, f[0], f[1], f[2], f[3], f[4], f[5], (n >= 13 ? f + 6 : NULL));
s_RegisterUndoForAtomAttrChange(self, s_AnisoSym, val, oval);
return val;
}
{"fract_x", &s_FractXSym, 0, s_AtomRef_GetFractionalX, s_AtomRef_SetFractionalX},
{"fract_y", &s_FractYSym, 0, s_AtomRef_GetFractionalY, s_AtomRef_SetFractionalY},
{"fract_z", &s_FractZSym, 0, s_AtomRef_GetFractionalZ, s_AtomRef_SetFractionalZ},
+ {"sigma", &s_SigmaSym, 0, s_AtomRef_GetSigma, s_AtomRef_SetSigma},
+ {"sigma_x", &s_SigmaXSym, 0, s_AtomRef_GetSigmaX, s_AtomRef_SetSigmaX},
+ {"sigma_y", &s_SigmaYSym, 0, s_AtomRef_GetSigmaY, s_AtomRef_SetSigmaY},
+ {"sigma_z", &s_SigmaZSym, 0, s_AtomRef_GetSigmaZ, s_AtomRef_SetSigmaZ},
{"v", &s_VSym, 0, s_AtomRef_GetV, s_AtomRef_SetV},
{"f", &s_FSym, 0, s_AtomRef_GetF, s_AtomRef_SetF},
{"occupancy", &s_OccupancySym, 0, s_AtomRef_GetOccupancy, s_AtomRef_SetOccupancy},
/*
* call-seq:
* cell = [a, b, c, alpha, beta, gamma]
- * set_cell([a, b, c, alpha, beta, gamma], flag = nil)
+ * set_cell([a, b, c, alpha, beta, gamma], flag = nil, sigma = nil)
*
* Set the unit cell parameters. If the cell value is nil, then clear the current cell.
This operation is undoable. If the second argument is given as non-nil, then
the coordinates are transformed so that the cartesian coordinates remain the same.
+ If the third argument is given, then standard deviations of the cell is also set.
*/
static VALUE
s_Molecule_SetCell(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE val, fval;
+ VALUE val, fval, sval;
int i, flag;
- double d[6];
+ double d[12];
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "11", &val, &fval);
+ rb_scan_args(argc, argv, "12", &val, &fval, &sval);
flag = RTEST(fval);
if (val == Qnil) {
MolActionCreateAndPerform(mol, gMolActionSetCell, 0, d, flag);
}
for (i = 0; i < 6; i++)
d[i] = NUM2DBL(rb_Float(Ruby_ObjectAtIndex(val, i)));
- MolActionCreateAndPerform(mol, gMolActionSetCell, 6, d, flag);
+ if (sval != Qnil) {
+ for (i = 6; i < 11; i++)
+ d[i] = NUM2DBL(rb_Float(Ruby_ObjectAtIndex(sval, i - 6)));
+ }
+ MolActionCreateAndPerform(mol, gMolActionSetCell, i, d, flag);
// MoleculeSetCell(mol, d[0], d[1], d[2], d[3], d[4], d[5]);
return val;
}
return @tokens.shift
end
def float_strip_rms(str)
- return Float(str.sub(/\(\d+\)/, ""))
+ str =~ /^(-?)(\d*)(\.(\d*))?(\((\d+)\))?$/
+ sgn, i, frac, rms = $1, $2, $4, $6
+ i = i.to_f
+ if frac
+ base = 0.1 ** frac.length
+ i = i + frac.to_f * base
+ else
+ base = 1.0
+ end
+ if rms
+ rms = rms.to_f * base
+ else
+ rms = 0.0
+ end
+ if sgn == "-"
+ i = -i
+ end
+ return i, rms
end
@tokens = []
self.remove(All)
if token =~ /^_cell/
val = getciftoken(fp)
if token == "_cell_length_a"
- cell[0] = float_strip_rms(val)
+ cell[0], cell[6] = float_strip_rms(val)
elsif token == "_cell_length_b"
- cell[1] = float_strip_rms(val)
+ cell[1], cell[7] = float_strip_rms(val)
elsif token == "_cell_length_c"
- cell[2] = float_strip_rms(val)
+ cell[2], cell[8] = float_strip_rms(val)
elsif token == "_cell_angle_alpha"
- cell[3] = float_strip_rms(val)
+ cell[3], cell[9] = float_strip_rms(val)
elsif token == "_cell_angle_beta"
- cell[4] = float_strip_rms(val)
+ cell[4], cell[10] = float_strip_rms(val)
elsif token == "_cell_angle_gamma"
- cell[5] = float_strip_rms(val)
+ cell[5], cell[11] = float_strip_rms(val)
end
- if cell.length == 6 && cell.all?
+ if cell.length == 12 && cell.all?
self.cell = cell
puts "Unit cell is set to #{cell.inspect}."
cell = []
occ = d[hlabel["_atom_site_occupancy"]]
calc = d[hlabel["_atom_site_calc_flag"]]
ap = self.add_atom(name, elem, elem)
- ap.fract_x = float_strip_rms(fx)
- ap.fract_y = float_strip_rms(fy)
- ap.fract_z = float_strip_rms(fz)
+ ap.fract_x, ap.sigma_x = float_strip_rms(fx)
+ ap.fract_y, ap.sigma_y = float_strip_rms(fy)
+ ap.fract_z, ap.sigma_z = float_strip_rms(fz)
if biso
- ap.temp_factor = float_strip_rms(biso)
+ ap.temp_factor, sig = float_strip_rms(biso)
elsif uiso
- ap.temp_factor = float_strip_rms(uiso) * 78.9568352087149 # 8*pi*pi
+ ap.temp_factor, sig = float_strip_rms(uiso) * 78.9568352087149 # 8*pi*pi
end
- ap.occupancy = float_strip_rms(occ)
+ ap.occupancy, sig = float_strip_rms(occ)
if calc == "c" || calc == "calc"
calculated_atoms.push(ap.index)
end
next if !ap
u11 = d[hlabel["_atom_site_aniso_U_11"]]
if u11
- u11 = float_strip_rms(u11)
- u22 = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
- u33 = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
- u12 = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
- u13 = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
- u23 = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
- ap.aniso = [u11, u22, u33, u12, u13, u23, 8]
+ usig = []
+ u11, usig[0] = float_strip_rms(u11)
+ u22, usig[1] = float_strip_rms(d[hlabel["_atom_site_aniso_U_22"]])
+ u33, usig[2] = float_strip_rms(d[hlabel["_atom_site_aniso_U_33"]])
+ u12, usig[3] = float_strip_rms(d[hlabel["_atom_site_aniso_U_12"]])
+ u13, usig[4] = float_strip_rms(d[hlabel["_atom_site_aniso_U_13"]])
+ u23, usig[5] = float_strip_rms(d[hlabel["_atom_site_aniso_U_23"]])
+ ap.aniso = [u11, u22, u33, u12, u13, u23, 8] + usig
c += 1
end
}
OSDN Git Service
