-
Notifications
You must be signed in to change notification settings - Fork 0
/
conformation.lisp
167 lines (150 loc) · 9.88 KB
/
conformation.lisp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
(in-package :topology)
(defclass conformation ()
((monomer-positions :initarg :monomer-positions :accessor monomer-positions)
(monomer-contexts :type hash-table :initarg :monomer-contexts :accessor monomer-contexts)
(oligomers :initarg :oligomers :accessor oligomers)
(aggregate :initarg :aggregate :accessor aggregate)
(energy-function :initarg :energy-function :accessor energy-function)
(ataggregate :initarg :ataggregate :accessor ataggregate)
(joint-tree :initarg :joint-tree :accessor joint-tree)))
(cando.serialize:make-class-save-load conformation)
(defgeneric foldamer-monomer-context (focus-monomer oligomer foldamer)
(:documentation "Return a monomer-context for a monomer in the oligomer using the foldamer.
Specialize the foldamer argument to provide methods"))
(defun make-conformation (oligomers &key monomer-order)
"Build a conformation for the oligomers."
(let* ((aggregate (chem:make-aggregate :all))
(monomer-positions (make-hash-table))
(oligomer-molecules (loop for oligomer in oligomers
for molecule = (topology:build-molecule oligomer aggregate monomer-positions)
collect (cons oligomer molecule)))
(monomer-contexts (make-hash-table))
(ataggregate (let ((atagg (make-instance 'ataggregate :aggregate aggregate)))
(resize-atmolecules atagg (length oligomers))
atagg))
(joint-tree (make-joint-tree))
(energy-function (chem:make-energy-function :matter aggregate))
)
(loop for oligomer-molecule in oligomer-molecules
for oligomer = (car oligomer-molecule)
for molecule = (cdr oligomer-molecule)
for oligomer-space = (oligomer-space oligomer)
for foldamer = (foldamer oligomer-space)
for molecule-index from 0
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 molecule molecule-index monomer-positions joint-tree (chem:atom-table energy-function))))
(put-atmolecule ataggregate atmolecule molecule-index))
finally (return (make-instance 'conformation
:monomer-positions monomer-positions
:monomer-contexts monomer-contexts
:oligomers (list oligomer)
:aggregate aggregate
:energy-function energy-function
:ataggregate ataggregate
:joint-tree joint-tree)))
))
(defun walk-atoms-joints (conformation callback)
(let ((aggregate (aggregate conformation))
(ataggregate (ataggregate conformation)))
(loop for molecule-index below (chem:content-size aggregate)
for molecule = (chem:content-at aggregate molecule-index)
for atmolecule = (aref (atmolecules ataggregate) molecule-index)
do (loop for residue-index below (chem:content-size molecule)
for residue = (chem:content-at molecule residue-index)
for atresidue = (aref (atresidues atmolecule) residue-index)
do (loop for atom-index below (chem:content-size residue)
for atom = (chem:content-at residue atom-index)
for joint = (aref (joints atresidue) atom-index)
do (funcall callback atom joint (list molecule-index residue-index atom-index)))))))
#+(or)
(defun copy-atom-positions-into-joints (conformation)
(walk-atoms-joints conformation
(lambda (atm jnt atomid)
(declare (ignore atomid))
(kin:set-position jnt (chem:get-position atm)))))
#+(or)(defun copy-joint-positions-into-atoms (conformation)
(walk-atoms-joints conformation
(lambda (atm jnt atomid)
(declare (ignore atomid))
(chem:set-position atm (kin:position coords jnt)))))
(defun update-joint-tree-internal-coordinates (conformation coordinates)
(let ((ataggregate (ataggregate conformation)))
(walk-ataggregate-joints ataggregate
(lambda (joint atom-id)
(declare (ignore atom-id))
(kin:update-internal-coord joint coordinates)))))
(defun build-all-atom-tree-external-coordinates (conformation coords)
(let ((joint (root (joint-tree conformation))))
(kin:update-xyz-coords joint coords)))
(defun zero-all-atom-tree-external-coordinates (conf)
"Set the external coordinates for each joint to the origin"
(let ((ataggregate (ataggregate conf)))
(walk-ataggregate-joints ataggregate
(lambda (joint atom-id)
(declare (ignore atom-id))
(kin:set-position joint (geom:vec 0.0 0.0 0.0))))))
(defun fill-internals-from-oligomer-shape (conf fragments oligomer-shape)
"Fill internal coordinates from the fragments"
(loop for oligomer in (oligomers conf)
do (loop with atagg = (ataggregate conf)
with atmol = (elt (atmolecules atagg) 0)
for monomer in (ordered-monomers oligomer)
for monomer-context = (gethash monomer (monomer-contexts conf))
for fragment-conformations = (let ((frag (gethash monomer-context (monomer-context-to-fragment-conformations fragments))))
(unless frag
(error "Could not find monomer-context ~a" monomer-context))
frag)
for monomer-index = (gethash monomer (monomer-positions conf))
for atres = (elt (atresidues atmol) monomer-index)
for monomer-shape = (let ((ms (gethash monomer (monomer-shape-map oligomer-shape))))
(unless ms (error "Could not get monomer-shape for monomer ~a" monomer))
ms)
for fragment-conformation-index = (let ((fi (fragment-conformation-index monomer-shape)))
(unless fi (break "Could not find fragment-conformation-index ~a ~a" fragment-conformation-index monomer-shape))
fi)
for fragment-internals = (let ((fragments (fragments fragment-conformations)))
(unless fragments
(error "fragments is NIL for context ~a" monomer-context))
(elt (fragments fragment-conformations) fragment-conformation-index))
do (loop for joint across (joints atres)
for internal across (internals fragment-internals)
do (fill-joint-internals joint internal))
)))
(defun search-conformations (oligomer-space fragment-conformations monomer-names sdf-filename &optional (number-struct 50) (number-conf 100))
(error "Implement search-conformations")
#+(or)(with-open-file (fout sdf-filename :direction :output)
(loop for struct-count below number-struct
for rand-sequence = (random (topology:number-of-sequences oligomer-space))
for _a = (format t "rand-sequence ~a~%" rand-sequence)
for oligomer = (topology:make-oligomer oligomer-space rand-sequence)
for conf = (topology:make-conformation oligomer)
do (topology::fill-internals-from-fragments conf fragment-conformations 0)
(loop for count below number-conf
do (topology::fill-internals-from-fragments-for-monomers-named conf fragment-conformations monomer-names)
do (topology:zero-all-atom-tree-external-coordinates conf)
do (topology:build-all-atom-tree-external-coordinates conf)
do (topology:copy-joint-positions-into-atoms conf)
do (sdf:write-sdf-stream (topology:aggregate conf) fout)))))
(defun fill-internals-from-fragments-for-monomers-named (conf fragments monomer-names)
"Fill internal coordinates from the fragments"
(loop for oligomer in (oligomers conf)
do (loop with atagg = (ataggregate conf)
with atmol = (elt (atmolecules atagg) 0)
for monomer across (monomers oligomer)
when (member (topology:current-stereoisomer-name monomer oligomer) monomer-names)
do (let* ((monomer-context (gethash monomer (monomer-contexts conf)))
(frags (gethash monomer-context (monomer-context-to-fragment-conformations fragments)))
(monomer-index (gethash monomer (monomer-positions conf)))
(atres (elt (atresidues atmol) monomer-index))
(fragment-conformations (gethash monomer-context (monomer-context-to-fragment-conformations fragments)))
(rand-limit (length (fragments fragment-conformations)))
(rand-index (random rand-limit))
(fragment-internals (elt (fragments fragment-conformations) rand-index)))
(loop for joint across (joints atres)
for internal in (internals fragment-internals)
do (fill-joint-internals joint internal))
))))