Switch to ENSURE_NOT_NAN

include/cando/chem/energyComponent.h

@@ -179,7 +179,7 @@ inline string	atomLabel(Atom_sp a)
 
 #define	ForceAcc(i,o,v) {\
       if ( hasForce ) {\
-        ASSERT_NOT_NAN(v); \
+        ENSURE_NOT_NAN(v); \
         force->setElement((i)+(o),(v)+force->getElement((i)+(o)));\
       }\
   }
@@ -189,14 +189,14 @@ inline string	atomLabel(Atom_sp a)
 //
 #define	OffDiagHessAcc(i1,o1,i2,o2,v) {\
     if ( hasHessian ) {\
-      ASSERT_NOT_NAN(v); \
+      ENSURE_NOT_NAN(v); \
       hessian->addToElement((i1)+(o1),(i2)+(o2),v);\
     }\
     if ( hasHdAndD ) {\
       auto v22 = v*dvec->element((i2)+(o2));\
       auto v11 = v*dvec->element((i1)+(o1));\
-      ASSERT_NOT_NAN(v22);\
-      ASSERT_NOT_NAN(v11);\
+      ENSURE_NOT_NAN(v22);\
+      ENSURE_NOT_NAN(v11);\
       hdvec->addToElement((i1)+(o1),v22); \
       hdvec->addToElement((i2)+(o2),v11); \
     }\
@@ -207,12 +207,12 @@ inline string	atomLabel(Atom_sp a)
 //
 #define	DiagHessAcc(i1,o1,i2,o2,v) {\
     if ( hasHessian ) {\
-      ASSERT_NOT_NAN(v); \
+      ENSURE_NOT_NAN(v); \
       hessian->addToElement((i1)+(o1),(i2)+(o2),v);\
     }\
     if ( hasHdAndD ) {\
       auto vd = v*dvec->element((i1)+(o1));\
-      ASSERT_NOT_NAN(vd); \
+      ENSURE_NOT_NAN(vd); \
       hdvec->addToElement((i1)+(o1),vd); \
     }\
   }

src/chem/cluster.cc

@@ -114,7 +114,7 @@ int   Kmeans_O::NearestCenter(core::SimpleVector_sp centers, Point p)
   for (int k = 0; k < _K; k++ )
   {
     dis = Distance<floatType>(p, gc::As<Point>((*centers)[k]) );
-    ASSERT_NOT_NAN(dis);
+    ENSURE_NOT_NAN(dis);
     if (dis < minDistance) {
       minDistance = dis;
       k_id = k;

src/chem/energyAngle.cc

@@ -402,12 +402,12 @@ bool		calcOffDiagonalHessian = true;
 #define	ANGLE_FORCE_ACCUMULATE(i,o,v) {}
 #undef	ANGLE_DIAGONAL_HESSIAN_ACCUMULATE
 #define	ANGLE_DIAGONAL_HESSIAN_ACCUMULATE(i1,o1,i2,o2,v) {\
-    ASSERT_NOT_NAN(v); \
+    ENSURE_NOT_NAN(v); \
     m->addToElement((i1)+(o1),(i2)+(o2),v);\
 }
 #undef	ANGLE_OFF_DIAGONAL_HESSIAN_ACCUMULATE
 #define	ANGLE_OFF_DIAGONAL_HESSIAN_ACCUMULATE(i1,o1,i2,o2,v) {\
-    ASSERT_NOT_NAN(v); \
+    ENSURE_NOT_NAN(v); \
     m->addToElement((i1)+(o1),(i2)+(o2),v);\
 }
 #define ANGLE_CALC_FORCE
@@ -483,7 +483,7 @@ double EnergyAngle_O::evaluateAllComponent( ScoringFunction_sp score,
 #undef ANGLE_SET_POSITION
 #define ANGLE_SET_POSITION(x,ii,of)	{x=pos->element(ii+of);}
 #undef ANGLE_ENERGY_ACCUMULATE
-#define ANGLE_ENERGY_ACCUMULATE(e) { ASSERT_NOT_NAN(e); termEnergy += (e); }
+#define ANGLE_ENERGY_ACCUMULATE(e) { ENSURE_NOT_NAN(e); termEnergy += (e); }
 #undef	ANGLE_FORCE_ACCUMULATE
 #undef	ANGLE_DIAGONAL_HESSIAN_ACCUMULATE
 #undef	ANGLE_OFF_DIAGONAL_HESSIAN_ACCUMULATE

src/chem/energyDihedral.cc

@@ -777,7 +777,7 @@ if (hasActiveAtomMask \
 #undef	DIHEDRAL_SET_POSITION
 #define	DIHEDRAL_SET_POSITION(x,ii,of)	{x=pos->element(ii+of);}
 #undef	DIHEDRAL_ENERGY_ACCUMULATE
-#define	DIHEDRAL_ENERGY_ACCUMULATE(e) {ASSERT_NOT_NAN(e); totalEnergy += (e); }
+#define	DIHEDRAL_ENERGY_ACCUMULATE(e) {ENSURE_NOT_NAN(e); totalEnergy += (e); }
 #undef	DIHEDRAL_FORCE_ACCUMULATE
 #undef	DIHEDRAL_DIAGONAL_HESSIAN_ACCUMULATE
 #undef	DIHEDRAL_OFF_DIAGONAL_HESSIAN_ACCUMULATE

src/chem/energyNonbond.cc

@@ -350,9 +350,9 @@ double template_evaluateUsingExcludedAtoms(EnergyNonbond_O* mthis,
 #undef	NONBOND_SET_POSITION
 #define	NONBOND_SET_POSITION(x,ii,of)	{x=pos->element(ii+of);}
 #undef	NONBOND_EEEL_ENERGY_ACCUMULATE
-#define	NONBOND_EEEL_ENERGY_ACCUMULATE(e) {ASSERT_NOT_NAN(e); energyElectrostatic +=(e);}
+#define	NONBOND_EEEL_ENERGY_ACCUMULATE(e) {ENSURE_NOT_NAN(e); energyElectrostatic +=(e);}
 #undef	NONBOND_EVDW_ENERGY_ACCUMULATE
-#define	NONBOND_EVDW_ENERGY_ACCUMULATE(e) {ASSERT_NOT_NAN(e); energyVdw+=(e);}
+#define	NONBOND_EVDW_ENERGY_ACCUMULATE(e) {ENSURE_NOT_NAN(e); energyVdw+=(e);}
 #undef	NONBOND_ENERGY_ACCUMULATE
 #define	NONBOND_ENERGY_ACCUMULATE(e) {};
 #undef	NONBOND_FORCE_ACCUMULATE
@@ -577,9 +577,9 @@ double template_evaluateUsingTerms(EnergyNonbond_O* mthis,
 #undef	NONBOND_SET_POSITION
 #define	NONBOND_SET_POSITION(x,ii,of)	{x=pos->element(ii+of);}
 #undef	NONBOND_EEEL_ENERGY_ACCUMULATE
-#define	NONBOND_EEEL_ENERGY_ACCUMULATE(e) {ASSERT_NOT_NAN(e); if (InteractionIs14) {energyElectrostatic14 +=(e);} else {energyElectrostatic +=(e);}}
+#define	NONBOND_EEEL_ENERGY_ACCUMULATE(e) {ENSURE_NOT_NAN(e); if (InteractionIs14) {energyElectrostatic14 +=(e);} else {energyElectrostatic +=(e);}}
 #undef	NONBOND_EVDW_ENERGY_ACCUMULATE
-#define	NONBOND_EVDW_ENERGY_ACCUMULATE(e) {ASSERT_NOT_NAN(e); if (InteractionIs14) {energyVdw14+=(e);} else {energyVdw+=(e);}}
+#define	NONBOND_EVDW_ENERGY_ACCUMULATE(e) {ENSURE_NOT_NAN(e); if (InteractionIs14) {energyVdw14+=(e);} else {energyVdw+=(e);}}
 
 #undef	NONBOND_ENERGY_ACCUMULATE
 #define	NONBOND_ENERGY_ACCUMULATE(e) {};

src/chem/energyStretch.cc

@@ -310,12 +310,12 @@ void	EnergyStretch_O::setupHessianPreconditioner(NVector_sp nvPosition,
 #define	STRETCH_FORCE_ACCUMULATE(i,o,v) {}
 #undef	STRETCH_DIAGONAL_HESSIAN_ACCUMULATE
 #define	STRETCH_DIAGONAL_HESSIAN_ACCUMULATE(i1,o1,i2,o2,v) {	\
-    ASSERT_NOT_NAN(v); \
+    ENSURE_NOT_NAN(v); \
     m->addToElement((i1)+(o1),(i2)+(o2),v);		\
   }
 #undef	STRETCH_OFF_DIAGONAL_HESSIAN_ACCUMULATE
 #define	STRETCH_OFF_DIAGONAL_HESSIAN_ACCUMULATE(i1,o1,i2,o2,v) {	\
-    ASSERT_NOT_NAN(v); \
+    ENSURE_NOT_NAN(v); \
     m->addToElement((i1)+(o1),(i2)+(o2),v);			\
   }
 #define STRETCH_CALC_FORCE
@@ -396,7 +396,7 @@ double EnergyStretch_O::evaluateAllComponent( ScoringFunction_sp score,
 #undef	STRETCH_SET_POSITION
 #define	STRETCH_SET_POSITION(x,ii,of)	{x = pos->getElement(ii+of);}
 #undef	STRETCH_ENERGY_ACCUMULATE
-#define	STRETCH_ENERGY_ACCUMULATE(e) {ASSERT_NOT_NAN(e);  totalEnergy += (e); }
+#define	STRETCH_ENERGY_ACCUMULATE(e) {ENSURE_NOT_NAN(e);  totalEnergy += (e); }
 #undef	STRETCH_FORCE_ACCUMULATE
 #undef	STRETCH_DIAGONAL_HESSIAN_ACCUMULATE
 #undef	STRETCH_OFF_DIAGONAL_HESSIAN_ACCUMULATE

src/chem/nVector.cc

@@ -436,9 +436,9 @@ CL_DEFUN vecreal chem__nvector_magnitude(NVector_sp vec) {
     val = val*val;
     sum2 += val;
   }
-  ASSERT_NOT_NAN(sum2);
+  ENSURE_NOT_NAN(sum2);
   vecreal res = sqrt(sum2);
-  ASSERT_NOT_NAN(res);
+  ENSURE_NOT_NAN(res);
   return res;
 }
 
@@ -655,7 +655,7 @@ CL_DEFUN vecreal chem__nvector_dot(NVector_sp veca, NVector_sp vecb) {
     val = vala*valb;
     sum2 += val;
   }
-  ASSERT_NOT_NAN(sum2);
+  ENSURE_NOT_NAN(sum2);
   return sqrt(sum2);
 }
 

src/geom/matrix.cc

@@ -38,6 +38,7 @@ at mailto:[email protected] if you would like a different license.
 
 #include <string.h>
 #include <cando/geom/matrix.h>
+#include <cando/geom/ovector3.h>
 #include <iostream>
 #include <clasp/core/foundation.h>
 #include <clasp/core/object.h>
@@ -1348,7 +1349,7 @@ void internalCoordinatesFromPointAndStub(const Vector3 &D,
   double dz = d.dotProduct(z);
   VEC_LOG("dx = {}  dy = {}  dz = {}\n", dx, dy, dz);
   phi = geom::geom__planeVectorAngle(dy, dz);
-  if (std::isnan(phi)) core::lisp_nan_error();
+  ENSURE_NOT_NAN(phi);
   VEC_LOG("  dy = {}   dz = {}\n", dy, dz);
   VEC_LOG("_Phi = {} deg\n", (phi / 0.0174533));
   Vector3 dox(1.0, 0.0, 0.0);
@@ -1374,7 +1375,7 @@ void internalCoordinatesFromPointAndStub(const Vector3 &D,
   VEC_LOG("eox = {}  eoy = {}\n", eox, eoy);
   //  double eoz = dop.dotProduct(doz); // Must be 0.0
   theta = - geom::geom__planeVectorAngle(eox, eoy);
-  if (std::isnan(theta)) core::lisp_nan_error();
+  ENSURE_NOT_NAN(theta);
   if (theta<0.0) {
     SIMPLE_WARN("Theta {} should never be less than zero", theta);
   }
@@ -1403,11 +1404,20 @@ void stubFromThreePoints(
   LOG(("a = %s") , a.asString());
   LOG(("b = %s") , b.asString());
   LOG(("c = %s") , c.asString());
+
   Vector3 e1(b-c);
-  e1 = e1.normalized();
+  if (e1.length()>0.0) {
+    e1 = e1.normalized();
+  } else {
+    CLASP_ERROR("When normalizing b-c it is zero length b=~s c=~s", OVector3_O::create(b), OVector3_O::create(c));
+  }
   LOG(("e1 = (b-c).normalized : %s") , e1.asString() );
   Vector3 e3(e1.crossProduct(a-c));
+  if (e3.length()>0.0) {
   e3 = e3.normalized();
+  } else {
+    CLASP_ERROR("When normalizing e3 it was zero length e1=~s a=~s c=~s", OVector3_O::create(e1), OVector3_O::create(a), OVector3_O::create(c));
+  }
   LOG(("e3 = (e1.crossProduct(a-c)).normalized : %s") , e3.asString() );
   Vector3 e2(e3.crossProduct(e1));
   LOG(("e2 = (e3.crossProduct(e1)): %s") , e2.asString() );

src/lisp/cando-widgets/show.lisp

@@ -534,7 +534,15 @@
     agg))
 
 (defmethod show-on-pane (pane-instance (object chem:aggregate) &rest rest &key &allow-other-keys)
-  (apply #'ngl-show-on-pane pane-instance object rest))
+  (let ((rev-molecules nil))
+    (chem:do-molecules (mol object)
+      (push mol rev-molecules))
+    (let ((molecules (nreverse rev-molecules)))
+      (loop for mol in molecules
+            for append = nil then t
+            with pane-instance = nil
+            do (setf pane-instance (apply #'show-on-pane pane-instance mol :append append rest))))))
+
 
 (defmethod show-on-pane (pane-instance (object chem:molecule) &rest rest &key &allow-other-keys)
   (apply #'ngl-show-on-pane pane-instance object rest))

src/lisp/topology/conformation.lisp

@@ -8,7 +8,7 @@
            :documentation "Apply this to the coordinates after the to-origin transform was applied")
    (from-origin :initarg :from-origin :accessor from-origin
                 :documentation "This is the second transform that takes the object after to-origin is applied"))
-  (:documentation "Represent the orientation of a molecule in a design calculation."))"
+  (:documentation "Represent the orientation of a molecule in a design calculation."))
 
 (defun make-orientation (&key (from-origin (geom:make-matrix-identity))
                            (adjust (geom:make-matrix-identity))
@@ -40,6 +40,7 @@
 
 (defmethod kin:orientation-transform ((orientation orientation))
   (let* ((from-origin (topology:from-origin orientation))
+         (adjust (adjust orientation))
          (to-origin (to-origin orientation))
          (transform0 (geom:m*m adjust to-origin))
          (transform (geom:m*m from-origin transform0)))

src/lisp/topology/jupyter.lisp

@@ -164,4 +164,4 @@
          (coords (topology:make-coordinates-for-assembler assembler)))
     (topology:fill-internals-from-oligomer-shape assembler object)
     (topology:build-external-coordinates assembler :coords coords)
-    (cando-widgets::show-on-pane pane-instance (topology:aggregate* assembler coords))))
+    (cando-widgets::show-on-pane pane-instance (cando:mol (topology:aggregate* assembler coords) 0))))

src/lisp/topology/oligomer.lisp

@@ -313,15 +313,16 @@ of monomers-to-residues."
   (let ((root-monomer (when (in-monomer-subset monomer-subset (root-monomer oligomer))))
         (monomer-out-couplings (make-hash-table))
         (monomers-to-topologys (make-hash-table))
-        (ring-couplings nil))
+        (ring-couplings nil)
+        (name (name (oligomer-space oligomer))))
     (declare (ignore root-monomer))
     (loop for index below (length (couplings oligomer))
           for coupling = (elt (couplings oligomer) index)
           if (typep coupling 'directional-coupling)
             do (push coupling (gethash (source-monomer coupling) monomer-out-couplings))
           else
             do (pushnew coupling ring-couplings)) ; Only add ring coupling when unique
-    (let* ((molecule (let ((mol (chem:make-molecule :mol)))
+    (let* ((molecule (let ((mol (chem:make-molecule name)))
                        (chem:setf-force-field-name mol (oligomer-force-field-name (foldamer (oligomer-space oligomer))))
                        mol)))
       (recursively-build-molecule oligomer

src/lisp/topology/packages.lisp

@@ -447,7 +447,8 @@
    #:make-rotamer-indexes
    #:make-indexes-into-rotamer-indexes
    #:rotamer-indexes
-   #:indexes-into-rotamer-indexes))
+   #:indexes-into-rotamer-indexes
+   #:ensure-permissible-rotamers-equal))
 
 (defpackage #:topology.dag
   (:use #:common-lisp)

src/lisp/topology/steppers.lisp

@@ -67,6 +67,15 @@
                                                                       (allowed-rotamer-indexes permissible-rotamer)))
                                 (coerce (permissible-rotamer-vector obj) 'list)))))
 
+(defun ensure-permissible-rotamers-equal (perm1 perm2)
+  (loop for p1 across (permissible-rotamer-vector perm1)
+        for p2 across (permissible-rotamer-vector perm2)
+        unless (equalp (allowed-rotamer-indexes p1)
+                       (allowed-rotamer-indexes p2))
+          do (error "perm1 ~s and~%perm2 ~s~%are not equalp" perm1 perm2))
+  )
+
+
 (defclass permissible-backbone-rotamers (permissible-rotamers)
   ())
 
@@ -371,6 +380,11 @@ Returns an instance of permissible-sidechain-rotamers."
   ((rotamer-indexes :initarg :rotamer-indexes :accessor rotamer-indexes))
   (:documentation "A class to index into a vector of rotamers"))
 
+(defmethod print-object ((object rotamer-indexes) stream)
+  (let ((*print-pretty* nil))
+    (print-unreadable-object (object stream :type t)
+      (format stream "~s" (rotamer-indexes object)))))
+
 (defun rotamer-indexes-length (rotamer-indexes)
   (length (rotamer-indexes rotamer-indexes)))
 
@@ -455,7 +469,8 @@ Returns an instance of permissible-sidechain-rotamers."
             for allowed-rotamer-indexes = (allowed-rotamer-indexes permissible-rotamer)
             for rotamer-index = (aref vec index)
             do (unless (position rotamer-index allowed-rotamer-indexes)
-                 (error "You are trying to write a rotamer-index ~a into locus ~a but it is not one of the permissible ones: ~s" rotamer-index locus allowed-rotamer-indexes))
+                 (let ((*print-pretty* nil))
+                   (error "You are trying to write a rotamer-index ~a into locus ~a but it is not one of the~%permissible ones: ~s" rotamer-index locus allowed-rotamer-indexes)))
             do (setf (topology:rotamer-index monomer-shape) rotamer-index))
       (ensure-oligomer-shape-is-consistent-with-permissible-rotamers oligomer-shape permissible-rotamers)
       oligomer-shape)))

src/lisp/topology/topology-classes.lisp

@@ -437,6 +437,7 @@ that is not avoid-out-coupling-plug-name.  Otherwise signal an error"
 
 (defclass oligomer-space ()
   ((foldamer :initarg :foldamer :accessor foldamer)
+   (name :initform :defos :initarg :name :reader name)
    (monomers :initform (make-array 16 :adjustable t :fill-pointer 0)
              :initarg :monomers :accessor monomers)
    (couplings :initform (make-array 16 :adjustable t :fill-pointer 0)
@@ -483,7 +484,7 @@ that is not avoid-out-coupling-plug-name.  Otherwise signal an error"
   (%number-of-sequences oligomer-space))
 
 
-(defun make-oligomer-space (foldamer tree &key (parts *parts*))
+(defun make-oligomer-space (foldamer tree &key (name :defos) (parts *parts*))
   "Make an oligomer-space from a description in the **tree**.
 The tree is a nested list of lists that look like 
 (component coupling component coupling component ... ).
@@ -499,7 +500,8 @@ Examples:
                  :default (cycle :first)))
 "
   (let* ((oligomer-space (make-instance 'oligomer-space
-                                        :foldamer foldamer))
+                                        :foldamer foldamer
+                                        :name name))
          (labels (make-hash-table)))
     (interpret-rooted-tree oligomer-space tree labels :parts parts)
     (setf (%number-of-sequences oligomer-space)