ORTEP silently fails with non positive-definite B. Draw isotropic atoms instead.

diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index dae904e..0fe2954 100755 (executable)
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -10284,6 +10284,7 @@ MoleculeSetAniso(Molecule *mp, int n1, int type, Double x11, Double x22, Double
        }
        MatrixSymDiagonalize(m1, val, axis);
        for (u = 0; u < 3; u++) {
+        anp->eigval[u] = val[u];
                if (val[u] < 0) {
                        fprintf(stderr, "Non-positive definite thermal parameters for atom %.4s\n", mp->atoms[n1].aname);
                        val[u] = 0.001;

diff --git a/MolLib/Molecule.h b/MolLib/Molecule.h
index 6afb0d5..716bc92 100755 (executable)
--- a/MolLib/Molecule.h
+++ b/MolLib/Molecule.h
@@ -46,6 +46,7 @@ typedef struct Aniso {
        char has_bsig;     /*  Has sigma values?  */
        Double  bsig[6];   /*  sigma values  */
        Mat33  pmat;      /*  A 3x3 matrix whose three column 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.  */
+    Double eigval[3]; /*  Eigenvalues of the B matrix; this should be all non-negative */
 } Aniso;
 
 /*  Symmetry operation  */

diff --git a/MolLib/Ruby_bind/ruby_bind.c b/MolLib/Ruby_bind/ruby_bind.c
index 6a9f2fc..57ac3ca 100644 (file)
--- a/MolLib/Ruby_bind/ruby_bind.c
+++ b/MolLib/Ruby_bind/ruby_bind.c
@@ -76,9 +76,9 @@ static VALUE
        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
@@ -3921,6 +3921,18 @@ static VALUE s_AtomRef_GetAniso(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);
@@ -4360,6 +4372,11 @@ static VALUE s_AtomRef_SetAniso(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;
@@ -4594,6 +4611,7 @@ static struct s_AtomAttrDef {
        {"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},

diff --git a/Scripts/crystal.rb b/Scripts/crystal.rb
index 62fb6f8..555f944 100755 (executable)
--- a/Scripts/crystal.rb
+++ b/Scripts/crystal.rb
@@ -206,6 +206,13 @@ def export_ortep(fp, attr = nil)
          }
        end
     an = ap.aniso
+    if an != nil
+      eigval = ap.aniso_eigenvalues
+      if eigval && (eigval[0] < 0.0 || eigval[1] < 0.0 || eigval[2] < 0.0)
+          #  Non positive-definite anisotropic factor: fallback to isotropic
+          an = nil
+      end
+    end
     if an != nil && rad < 0.0
       fp.printf " %8.5f%9.6f%9.6f%9.6f%9.6f%9.6f%9d\n", an[0], an[1], an[2], an[3], an[4], an[5], 0
     else
@@ -1531,8 +1538,11 @@ def cmd_show_ortep
          Dir.chdir(cwd)
          if pid != 0
            msg = "ORTEP execution in #{tmp} failed with status #{pid}."
-           message_box(msg, "ORTEP Failed", :ok, :warning)
-         else
+        message_box(msg, "ORTEP Failed", :ok, :warning)
+         elsif !FileTest.exist?(tmp + "/TEP001.ps")
+        msg = "ORTEP execution in #{tmp} failed with unknown error."
+        message_box(msg, "ORTEP Failed", :ok, :warning)
+      else
            open(tmp + "/TEP001.ps", "r") { |fp|
                  tepdata.clear
                  tepbounds = [100000, 100000, -100000, -100000]