-
Notifications
You must be signed in to change notification settings - Fork 17
/
topodihedrals.tcl
235 lines (199 loc) · 8.24 KB
/
topodihedrals.tcl
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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
#!/usr/bin/tclsh
# This file is part of TopoTools, a VMD package to simplify
# manipulating bonds other topology related properties.
#
# Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020 by Axel Kohlmeyer <[email protected]>
# $Id: topodihedrals.tcl,v 1.12 2020/05/29 19:47:40 johns Exp $
# return info about dihedrals
# we list and count only dihedrals that are entirely within the selection.
proc ::TopoTools::dihedralinfo {infotype sel {flag none}} {
set numdihedrals 0
array set dihedraltypes {}
set atomindex [$sel list]
set dihedrallist {}
foreach dihedral [join [molinfo [$sel molid] get dihedrals]] {
lassign $dihedral t a b c d
if {([lsearch -sorted -integer $atomindex $a] >= 0) \
&& ([lsearch -sorted -integer $atomindex $b] >= 0) \
&& ([lsearch -sorted -integer $atomindex $c] >= 0) \
&& ([lsearch -sorted -integer $atomindex $d] >= 0) } {
set dihedraltypes($t) 1
incr numdihedrals
lappend dihedrallist $dihedral
}
}
switch $infotype {
numdihedrals { return $numdihedrals }
numdihedraltypes { return [array size dihedraltypes] }
dihedraltypenames { return [lsort -ascii [array names dihedraltypes]] }
getdihedrallist { return $dihedrallist }
default { return "bug! shoot the programmer?"}
}
}
# delete all contained dihedrals of the selection.
proc ::TopoTools::cleardihedrals {sel} {
set mol [$sel molid]
set atomindex [$sel list]
set dihedrallist {}
foreach dihedral [join [molinfo $mol get dihedrals]] {
lassign $dihedral t a b c d
if {([lsearch -sorted -integer $atomindex $a] < 0) \
|| ([lsearch -sorted -integer $atomindex $b] < 0) \
|| ([lsearch -sorted -integer $atomindex $c] < 0) \
|| ([lsearch -sorted -integer $atomindex $d] < 0) } {
lappend dihedrallist $dihedral
}
}
molinfo $mol set dihedrals [list $dihedrallist]
}
# reset dihedrals to data in dihedrallist
proc ::TopoTools::setdihedrallist {sel dihedrallist} {
set mol [$sel molid]
set atomindex [$sel list]
set newdihedrallist {}
# set defaults
set t unknown; set a -1; set b -1; set c -1; set d -1
# preserve all dihedrals definitions that are not fully contained in $sel
foreach dihedral [join [molinfo $mol get dihedrals]] {
lassign $dihedral t a b c d
if {([lsearch -sorted -integer $atomindex $a] < 0) \
|| ([lsearch -sorted -integer $atomindex $b] < 0) \
|| ([lsearch -sorted -integer $atomindex $c] < 0) \
|| ([lsearch -sorted -integer $atomindex $d] < 0) } {
lappend newdihedrallist $dihedral
}
}
# append new ones, but only those contained in $sel
foreach dihedral $dihedrallist {
lassign $dihedral t a b c d
if {([lsearch -sorted -integer $atomindex $a] >= 0) \
&& ([lsearch -sorted -integer $atomindex $b] >= 0) \
&& ([lsearch -sorted -integer $atomindex $c] >= 0) \
&& ([lsearch -sorted -integer $atomindex $d] >= 0) } {
lappend newdihedrallist $dihedral
}
}
molinfo $mol set dihedrals [list $newdihedrallist]
}
# reset dihedrals to data in dihedrallist
proc ::TopoTools::retypedihedrals {sel} {
set mol [$sel molid]
set dihedrallist [dihedralinfo getdihedrallist $sel]
set atomtypes [$sel get type]
set atomindex [$sel list]
set newdihedrallist {}
foreach dihedral $dihedrallist {
lassign $dihedral type i1 i2 i3 i4
set idx [lsearch -sorted -integer $atomindex $i1]
set a [lindex $atomtypes $idx]
set idx [lsearch -sorted -integer $atomindex $i2]
set b [lindex $atomtypes $idx]
set idx [lsearch -sorted -integer $atomindex $i3]
set c [lindex $atomtypes $idx]
set idx [lsearch -sorted -integer $atomindex $i4]
set d [lindex $atomtypes $idx]
if { ([string compare $b $c] > 0) \
|| ( [string equal $b $c] && [string compare $a $d] > 0 ) } {
set t $a; set a $d; set d $t
set t $b; set b $c; set c $t
set t $i1; set i1 $i4; set i4 $t
set t $i2; set i2 $i3; set i3 $t
}
set type [join [list $a $b $c $d] "-"]
lappend newdihedrallist [list $type $i1 $i2 $i3 $i4]
}
setdihedrallist $sel $newdihedrallist
}
# reset dihedrals to definitions derived from bonds.
# this includes retyping of the dihedrals.
proc ::TopoTools::guessdihedrals {sel} {
set mol [$sel molid]
set atomtypes [$sel get type]
set atomindex [$sel list]
set newdihedrallist {}
set bondlist [bondinfo getbondlist $sel]
set bonddata [$sel getbonds]
# preserve all dihedrals definitions that are not fully contained in $sel
foreach dihedral [join [molinfo $mol get dihedrals]] {
lassign $dihedral t a b c d
if {([lsearch -sorted -integer $atomindex $a] < 0) \
|| ([lsearch -sorted -integer $atomindex $b] < 0) \
|| ([lsearch -sorted -integer $atomindex $c] < 0) \
|| ([lsearch -sorted -integer $atomindex $d] < 0) } {
lappend newdihedrallist $dihedral
}
}
# a topological dihedral is defined by a bond and atoms
# bound to it that are not the bond itself
foreach bond $bondlist {
lassign $bond b1 b2
set b1idx [lsearch -sorted -integer $atomindex $b1]
set b1typ [lindex $atomtypes $b1idx]
set b2idx [lsearch -sorted -integer $atomindex $b2]
set b2typ [lindex $atomtypes $b2idx]
foreach o1 [lindex $bonddata $b1idx] {
foreach o2 [lindex $bonddata $b2idx] {
if {($o1 == $b1) || ($o2 == $b1) || ($o1 == $b2) || ($o2 == $b2)} {
continue
}
set o1idx [lsearch -sorted -integer $atomindex $o1]
set o1typ [lindex $atomtypes $o1idx]
set o2idx [lsearch -sorted -integer $atomindex $o2]
set o2typ [lindex $atomtypes $o2idx]
if { $o1 != $o2 } {
if { ([string compare $b1typ $b2typ] > 0) \
|| ( [string equal $b1typ $b2typ]
&& [string compare $o1typ $o2typ] > 0 ) } {
set type [join [list $o2typ $b2typ $b1typ $o1typ] "-"]
lappend newdihedrallist [list $type $o2 $b2 $b1 $o1]
} else {
set type [join [list $o1typ $b1typ $b2typ $o2typ] "-"]
lappend newdihedrallist [list $type $o1 $b1 $b2 $o2]
}
}
}
}
}
setdihedrallist $sel $newdihedrallist
}
# define a new dihedral or change an existing one.
proc ::TopoTools::adddihedral {mol id1 id2 id3 id4 {type unknown}} {
if {[catch {atomselect $mol "index $id1 $id2 $id3 $id4"} sel]} {
vmdcon -err "topology adddihedral: Invalid atom indices: $sel"
return
}
# canonicalize indices
if {$id2 > $id3} {
set t $id2 ; set id2 $id3 ; set id3 $t
set t $id1 ; set id1 $id4 ; set id4 $t
}
set dihedrals [join [molinfo $mol get dihedrals]]
lappend dihedrals [list $type $id1 $id2 $id3 $id4]
molinfo $mol set dihedrals [list $dihedrals]
# this is not (yet) required
$sel delete
return
}
# delete a dihedral.
proc ::TopoTools::deldihedral {mol id1 id2 id3 id4 {type unknown}} {
if {[catch {atomselect $mol "index $id1 $id2 $id3 $id4"} sel]} {
vmdcon -err "topology deldihedral: Invalid atom indices: $sel"
return
}
# canonicalize indices
if {$id2 > $id3} {
set t $id2 ; set id2 $id3 ; set id3 $t
set t $id1 ; set id1 $id4 ; set id4 $t
}
set newdihedrallist {}
foreach dihedral [join [molinfo $mol get dihedrals]] {
lassign $dihedral t a b c d
if { ($a != $id1) || ($b != $id2) || ($c != $id3) || ($d != $id4) } {
lappend newdihedrallist $dihedral
}
}
molinfo $mol set dihedrals [list $newdihedrallist]
# this is not (yet) required
$sel delete
return
}