Skip to content

Commit

Permalink
Lots of topology changes
Browse files Browse the repository at this point in the history
  • Loading branch information
meister committed Nov 4, 2024
1 parent 13f2e81 commit 3fd3b2d
Show file tree
Hide file tree
Showing 16 changed files with 456 additions and 213 deletions.
1 change: 1 addition & 0 deletions include/cando/geom/omatrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ namespace geom {
public:
//! Create a 4x4 matrix uninitialized (identity=false) or identity if true
static OMatrix_sp make(bool identity);
CL_DEF_CLASS_METHOD static OMatrix_sp m4(core::List_sp args);
public:
string __repr__() const;

Expand Down
1 change: 1 addition & 0 deletions src/chem/matter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1192,6 +1192,7 @@ geom::Render_sp Matter_O::rendered(core::List_sp kopts)


CL_LISPIFY_NAME("allAtomsAsList");
CL_LAMBDA((matter matter) &optional allow-virtual-atoms);
CL_DEFMETHOD core::List_sp Matter_O::allAtomsAsList(bool allowVirtualAtoms ) const
{
core::List_sp result = nil<core::T_O>();
Expand Down
6 changes: 5 additions & 1 deletion src/chem/minimizer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -887,7 +887,11 @@ bool Minimizer_O::_displayIntermediateMessage(NVector_sp pos,
sout << fmt::format(" {:7}" , elapsed.count());
}
sout << fmt::format(" {:5}" , this->_Iteration);
sout << fmt::format(" {:9.4f}" , log(step));
if (step==0.0) {
sout << fmt::format(" {:9s}" , "nan" );
} else {
sout << fmt::format(" {:9.4f}" , log(step));
}
if ( steepestDescent )
{
sout << "StDesc";
Expand Down
9 changes: 9 additions & 0 deletions src/chem/nVector.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -190,6 +190,8 @@ double rmsMagnitudeWithActiveAtomMask(NVector_sp me, core::T_sp activeAtomMask)
}




double angleWithVector(NVector_sp me, NVector_sp other)
{
#define VERY_SMALL 1.0e-6
Expand Down Expand Up @@ -338,6 +340,13 @@ double rmsDistanceFromWithActiveAtomMask(NVector_sp u, NVector_sp v, core::T_sp
SIMPLE_ERROR("activeAtomMask must be a simple-bit-vector or NIL");
}

CL_DOCSTRING(R"doc(Calculate the root mean square difference between FIRST and SECOND.
If ACTIVE-MASK is defined then use it to select coordinates.)doc");
CL_LAMBDA(first second &optional active-mask);
CL_DEFUN double chem__nvector_rmsd(NVector_sp first, NVector_sp second, core::T_sp activeMask) {
return rmsDistanceFromWithActiveAtomMask(first,second,activeMask);
}


NVector_sp copy_nvector(NVector_sp source, size_t start, core::T_sp end )
{
Expand Down
22 changes: 22 additions & 0 deletions src/geom/omatrix.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,28 @@ CL_DEFUN OMatrix_sp OMatrix_O::make(bool ident)
return om;
};

OMatrix_sp OMatrix_O::m4(core::List_sp args)
{
if (core::cl__length(args)>16) {
SIMPLE_ERROR("The list of arguments {} is too long to fit in a 16 element homogeneous matrix", _rep_(args));
}
int index = 0;
auto om = gctools::GC<OMatrix_O>::allocate(true);
for ( auto cur : args ) {
core::T_sp obj = CONS_CAR(cur);
if (gc::IsA<core::Number_sp>(obj)) {
core::Number_sp num = gc::As_unsafe<core::Number_sp>(obj);
double val = core::clasp_to_double(num);
om->_Value[index] = val;
core::clasp_write_string(fmt::format("{} ", val ) );
} else {
TYPE_ERROR(obj,cl::_sym_Number_O);
}
index = index + 1;
}
return om;
};

string OMatrix_O::__repr__() const
{
stringstream ss;
Expand Down
5 changes: 5 additions & 0 deletions src/lisp/cando/convenience.lisp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,11 @@
(chem:content-at object 0)
(error "If no index is provided to mol there can only be one molecule and there are ~a molecules" (chem:content-size object))))))

(defmethod mol ((object chem:molecule) &optional id)
(if id
(error "When you pass a molecule to mol id must be nil - instead it is ~s" id)
object))

(defmethod res ((object chem:molecule) (id integer))
(chem:content-at object id))

Expand Down
2 changes: 1 addition & 1 deletion src/lisp/cando/modelling.lisp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@
"Use answer from here to get a random vector on a sphere.
https://math.stackexchange.com/questions/44689/how-to-find-a-random-axis-or-unit-vector-in-3d
"
(let* ((theta (random (* 2.0 (coerce PI 'double-flaot))))
(let* ((theta (random (* 2.0 (coerce PI 'double-float))))
(z (- (random 2.0) 1.0))
(sqrt-1-z2 (sqrt (- 1.0 (* z z))))
(x (* sqrt-1-z2 (cos theta)))
Expand Down
38 changes: 22 additions & 16 deletions src/lisp/cando/molecules.lisp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,18 +160,19 @@ Example: (set-stereoisomer-mapping *agg* '((:C1 :R) (:C2 :S))"
(chem:enable-print-intermediate-results minimizer)
(format t "enable-print-intermediate-results ~%"))
(chem:disable-print-intermediate-results minimizer))
(restart-case
(handler-bind
((chem:minimizer-error (lambda (err)
(warn "In minimize-no-fail - the minimizer reported: ~a" err)
(invoke-restart 'cando:skip-rest-of-minimization err))))
(with-handle-linear-angles-dihedrals (:max-times 3 :verbose verbose)
(chem:minimize minimizer active-atoms-mask)))
;; skip-rest-of-minimization can also be triggered by the user from the debugger
(skip-rest-of-minimization (err)
:report "Skip the rest of the current minimization - continue processing"
(chem:write-intermediate-results-to-energy-function minimizer)
(when resignal-error (error err)))))
(restart-case
(handler-bind
((chem:minimizer-error (lambda (err)
(warn "In minimize-no-fail - the minimizer reported: ~a" err)
(invoke-restart 'cando:skip-rest-of-minimization err))))
(with-handle-linear-angles-dihedrals (:max-times 3 :verbose verbose)
(ext:with-float-traps-masked (:underflow :overflow :invalid :inexact :divide-by-zero)
(chem:minimize minimizer active-atoms-mask))))
;; skip-rest-of-minimization can also be triggered by the user from the debugger
(skip-rest-of-minimization (err)
:report "Skip the rest of the current minimization - continue processing"
(chem:write-intermediate-results-to-energy-function minimizer)
(when resignal-error (error err)))))

(defmacro with-ignore-minimizer-errors (&body body)
`(handler-bind
Expand Down Expand Up @@ -214,7 +215,8 @@ Example: (set-stereoisomer-mapping *agg* '((:C1 :R) (:C2 :S))"
((ext:unix-signal-received (lambda (c)
(declare (ignore c))
(print "Done") (invoke-restart 'skip-rest))))
(chem:minimize minimizer))
(ext:with-float-traps-masked (:underflow :overflow :invalid :inexact :divide-by-zero)
(chem:minimize minimizer)))
(skip-rest ()
:report "Skip rest of minimization"
(chem:write-intermediate-results-to-energy-function minimizer)
Expand Down Expand Up @@ -736,7 +738,8 @@ Disabling happens before enabling - so you can disable all with T and then selec
(error "Exceeded max number of tries to build good geometry~%Bad geometry:~%~a" bad-geom)))


(defun save-mol2 (matter pathname &key use-sybyl-types)
(defun save-mol2 (matter pathname &key (use-sybyl-types t))
"Save MATTER as a mol2 to the file PATHNAME. USE-SYBYL-TYPES of nil requires atom types be set and then they will be used."
(let ((npn (merge-pathnames pathname *default-pathname-defaults*)))
(format t "Saving matter to ~a~%" npn)
(chem:save-mol2 matter npn use-sybyl-types)))
Expand All @@ -748,8 +751,11 @@ Disabling happens before enabling - so you can disable all with T and then selec
(chem:load-mol2 npn)))


(defun simple-vector-coordinate-for-atomspec (agg atomspec)
(let ((atoms (chem:atoms-with-chimera-specifications agg atomspec)))
(defun simple-vector-coordinate-for-atomspec (agg &optional atomspec)
"Return coordinates of AGG. If ATOMSPEC is not nil then use it to downselect atoms."
(let ((atoms (if atomspec
(chem:atoms-with-chimera-specifications agg atomspec)
(chem:matter/all-atoms-as-list agg))))
(chem:make-simple-vector-coordinate-from-atom-list atoms)))

(defun superpose-one (&key fixed-atoms moveable-matter moveable-atoms)
Expand Down
105 changes: 63 additions & 42 deletions src/lisp/topology/assembler.lisp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,8 +223,8 @@ Specialize the foldamer argument to provide methods"))
(t (error "Illegal value for monomer-subset ~s - must be NIL or a hash-table"))))

(defun make-assembler (oligomer-shapes &key orientations monomer-subset tune-energy-function (keep-interaction t))
"Build a assembler for the oligomers.
tune-energy-function - A function that takes the energy-function and an assembler and modifies the energy-function."
"Build a assembler for the OLIGOMER-SHAPES.
TUNE-ENERGY-FUNCTION - A function that takes the energy-function and an assembler and modifies the energy-function."
(cond
((not (or (= (length oligomer-shapes) 1) orientations))
(error "You must provide orientations when there is more than one oligomer-shape"))
Expand All @@ -236,46 +236,48 @@ tune-energy-function - A function that takes the energy-function and an assemble
(error "You must provide a list of oligomer-shapes"))
(let* ((aggregate (chem:make-aggregate :all))
(monomer-positions (make-hash-table))
(oligomers (mapcar (lambda (oligomer-shape)
(oligomer oligomer-shape))
oligomer-shapes))
(foldamers (mapcar (lambda (oligomer)
(let* ((oligomer-space (oligomer-space oligomer))
(foldamers (mapcar (lambda (oligomer-shape)
(let* ((oligomer (oligomer oligomer-shape))
(oligomer-space (oligomer-space oligomer))
(foldamer (foldamer oligomer-space)))
foldamer))
oligomers))
oligomer-shapes))
(monomers-to-residues (make-hash-table))
(oligomer-molecules (loop for oligomer in oligomers
for oligomer-space = (oligomer-space oligomer)
for foldamer = (foldamer oligomer-space)
for molecule-index from 0
for molecule = (topology:build-molecule oligomer
:monomer-subset monomer-subset
:aggregate aggregate
:molecule-index molecule-index
:monomers-to-residues monomers-to-residues
:monomer-positions-accumulator monomer-positions)
do (chem:setf-force-field-name molecule (oligomer-force-field-name foldamer))
collect (cons oligomer molecule))))
(oligomer-shapes-molecules (loop for oligomer-shape in oligomer-shapes
for oligomer = (oligomer oligomer-shape)
for oligomer-space = (oligomer-space oligomer)
for foldamer = (foldamer oligomer-space)
for molecule-index from 0
for molecule = (topology:build-molecule oligomer
:monomer-subset monomer-subset
:aggregate aggregate
:molecule-index molecule-index
:monomers-to-residues monomers-to-residues
:monomer-positions-accumulator monomer-positions)
do (chem:setf-force-field-name molecule (oligomer-force-field-name foldamer))
collect (cons oligomer-shape molecule))))
(let* ((monomer-contexts (make-hash-table))
(ataggregate (let ((atagg (make-instance 'ataggregate :aggregate aggregate)))
(resize-atmolecules atagg (length oligomers))
(resize-atmolecules atagg (length oligomer-shapes))
atagg))
(joint-tree (make-joint-tree))
(adjustments (make-instance 'adjustments))
(energy-function (chem:make-energy-function :keep-interaction keep-interaction :matter aggregate))
(assembler (loop for oligomer-molecule in oligomer-molecules
for oligomer = (car oligomer-molecule)
for molecule = (cdr oligomer-molecule)
(assembler (loop for oligomer-shape-molecule in oligomer-shapes-molecules
for oligomer-shape = (car oligomer-shape-molecule)
for oligomer = (oligomer oligomer-shape)
for molecule = (cdr oligomer-shape-molecule)
for oligomer-space = (oligomer-space oligomer)
for foldamer = (foldamer oligomer-space)
for molecule-index from 0
for monomer-to-monomer-shape-map = (monomer-shape-map oligomer-shape)
do (loop for monomer across (monomers oligomer-space)
for monomer-context = (foldamer-monomer-context monomer oligomer foldamer)
do (setf (gethash monomer monomer-contexts) monomer-context))
;; This is where I would invoke Conformation_O::buildMoleculeUsingOligomer
;; Use the monomers-to-topologys
do (let* ((atmolecule (build-atmolecule-using-oligomer oligomer
monomer-to-monomer-shape-map
monomer-subset
molecule
molecule-index
Expand Down Expand Up @@ -351,15 +353,15 @@ tune-energy-function - A function that takes the energy-function and an assemble
do (setf (gethash monomer monomer-contexts) monomer-context)))
;; This is where I would invoke Conformation_O::buildMoleculeUsingOligomer
;; Use the monomers-to-topologys
do (let* ((atmolecule (build-atmolecule-using-oligomer
oligomer
nil
molecule
molecule-index
monomer-positions
joint-tree
(chem:atom-table energy-function)
nil)))
do (let* ((atmolecule (build-atmolecule-using-oligomer oligomer
nil ; no monomer-to-monomer-shape-map
nil
molecule
molecule-index
monomer-positions
joint-tree
(chem:atom-table energy-function)
nil)))
(put-atmolecule ataggregate atmolecule molecule-index))
finally (return (make-instance 'training-assembler
:monomer-positions monomer-positions
Expand All @@ -373,8 +375,8 @@ tune-energy-function - A function that takes the energy-function and an assemble

(defun walk-atoms-joints (assembler molecule-index callback)
(let* ((aggregate (aggregate assembler))
(ataggregate (ataggregate assembler))
(molecule (chem:content-at aggregate molecule-index))
(ataggregate (ataggregate assembler))
(molecule (chem:content-at aggregate molecule-index))
(atmolecule (aref (atmolecules ataggregate) molecule-index)))
(loop for residue-index below (chem:content-size molecule)
for residue = (chem:content-at molecule residue-index)
Expand Down Expand Up @@ -460,12 +462,20 @@ tune-energy-function - A function that takes the energy-function and an assemble
(coords (topology:make-coordinates-for-assembler assembler)))
(if oligomer-shape
(progn
(unless orientationp
(when (and oligomer-shape (not orientationp))
(error "You must provide orientation when you provide oligomer-shape"))
(build-atom-tree-external-coordinates* assembler coords oligomer-shape orientation)
(adjust-atom-tree-external-coordinates assembler coords oligomer-shape))
(loop for oligomer-shape in (oligomer-shapes assembler)
do (build-external-coordinates assembler :oligomer-shape oligomer-shape :coords coords))))
do (build-external-coordinates assembler :oligomer-shape oligomer-shape
:orientation oligomer-shape
:coords coords))))

(defun update-externals (assembler &key oligomer-shape
(orientation :identity orientationp)
(coords (topology:make-coordinates-for-assembler assembler)))
(build-external-coordinates assembler :oligomer-shape oligomer-shape :orientation orientation :coords coords))


(defun build-atom-tree-external-coordinates* (assembler coords oligomer-shape maybe-orientation)
(let* ((orientation (orientation maybe-orientation assembler))
Expand Down Expand Up @@ -684,11 +694,7 @@ Some specialized methods will need coordinates for the assembler"))
monomers-to-residues))

(defmethod analyze-assembler ((assembler assembler) &optional name)
(loop for oligomer-shape in (oligomer-shapes assembler)
for oligomer-shape-index from 0
do (analyze-oligomer-shape oligomer-shape))
(let (
(backbone-internals-count (make-array (length (oligomer-shapes assembler)) :initial-element 0))
(let ((backbone-internals-count (make-array (length (oligomer-shapes assembler)) :initial-element 0))
(backbone-internals-defined (make-array (length (oligomer-shapes assembler)) :initial-element 0))
(sidechain-internals-count (make-array (length (oligomer-shapes assembler)) :initial-element 0))
(sidechain-internals-defined (make-array (length (oligomer-shapes assembler)) :initial-element 0))
Expand Down Expand Up @@ -732,3 +738,18 @@ Some specialized methods will need coordinates for the assembler"))
(format t "atres: ~s~%" atres)
(loop for joint in undefs
do (format t " joint: ~s undefined~%" joint))))))

(defun assembler-dump-internals (assembler)
(loop for oligomer-shapes in (oligomer-shapes assembler)
for index from 0
do (walk-atoms-joints assembler index
(lambda (atom joint atomid)
(format t "~s ~s~%" atom joint)))))


(defun assembler-dump-monomer-shapes (assembler)
(loop for oligomer-shape in (oligomer-shapes assembler)
for index from 0
do (format t "oligomer-shape ~s~%" oligomer-shape)
do (loop for monomer-shape across (topology:monomer-shape-vector oligomer-shape)
do (format t "monomer-shape ~s~%" monomer-shape))))
Loading

0 comments on commit 3fd3b2d

Please sign in to comment.