Fixed bug 417.
[gromacs/rigid-bodies.git] / scripts / xplor2gmx.pl
blob7769e030445521e43b7d618a83b84be84665b97c
1 #!/usr/bin/perl -w
4 # This script reads an XPLOR input file with distance restraint data
5 # as sometimes is found in the pdb database (http://www.pdb.org).
6 # From this input file dihedral restrints should be removed, such that
7 # only distance restraints are left. The script can handle ambiguous restraints.
8 # It converts the distance restraints to GROMACS format.
9 # A typical command line would be
10 # ./xplor2gmx.pl 33 conf.pdb < restraints.dat > disre.itp
11 # You can turn off debugging info, but please do verify your output is correct!
13 # David van der Spoel (spoel@gromacs.org), July 2002
16 # Turn debugging off (0), or on ( > 0).
17 $debug = 1;
18 # Turn atom name translation on and off
19 $trans = 1;
21 $res0 = shift;# || die "I need the residue offset\n";
22 $pdb = shift || die "I need the name of the pdb file with correct atom numbers\n";
23 $core = shift || "core.ndx";
24 $tbl = "/home/spoel/gmxdev/share/top/atom_nom.tbl";
26 printf "[ distance_restraints ]\n";
27 printf "; Read an xplor distance restraint file, and output GROMACS distance restraints.\n";
28 printf "; This also needs a pdb file with correct GROMACS atom numbering.\n";
29 printf "; I used $pdb for the atom numbers\n";
30 printf "; This also needs the offset in residues.\n";
31 printf "; I used $res0 for the residue offset\n";
33 # Read the pdb file
34 # If things go wrong, check whether your pdb file is read correctly.
35 $natom = 0;
36 $nres = 0;
37 open(PDB,$pdb) || die "Can not open file $pdb\n";
38 while ($line = <PDB>) {
39 if (index($line,"ATOM") >= 0) {
40 @tmp = split(' ',$line);
41 $aname[$natom] = $tmp[2];
42 $resnr[$natom] = $tmp[4];
43 if ($resnr[$natom] > $nres) {
44 $rname[$nres] = $tmp[3];
45 $nres++;
47 $natom++;
50 close PDB;
51 printf "; I found $natom atoms in the pdb file $pdb\n";
52 printf "; I found $nres residues in the pdb file $pdb\n";
53 if ($debug > 1) {
54 for ($i = 0; ($i < $natom); $i ++) {
55 printf("; %5d %5s %5s %5d\n",$i+1,$aname[$i],
56 $rname[$resnr[$i]-1],$resnr[$i]);
61 # Read the name translation table.
62 # Source for this comes from: http://www.bmrb.wisc.edu/ref_info/atom_nom.tbl
63 # It was adapted slightly for GROMACS use, but only as regards formatting,
64 # not for content
66 open(TBL,$tbl) || die "Can not open atom-name translation table $tbl\n";
67 $ntbl=0;
68 while ($line = <TBL>) {
69 @ttt = split('#',$line);
70 @tmp = split(' ',$ttt[0]);
71 if ($#tmp == 2) {
72 # New table entry
73 $tblres[$ntbl] = $tmp[0];
74 $tblxplor[$ntbl] = $tmp[1];
75 $tblgmx[$ntbl] = $tmp[2];
76 $ntbl++;
79 close TBL;
80 printf "; Read $ntbl entries from $tbl\n";
82 @templates = (
83 [ "HA#", "HA1", "HA2" ],
84 [ "HA*", "HA1", "HA2" ],
85 [ "HB#", "HB", "HB1", "HB2" ],
86 [ "HB*", "HB", "HB1", "HB2" ],
87 [ "HG#", "HG", "HG1", "HG2", "HG11", "HG12", "HG13", "HG21", "HG22", "HG23" ],
88 [ "HG*", "HG", "HG1", "HG2", "HG11", "HG12", "HG13", "HG21", "HG22", "HG23" ],
89 [ "HG1#", "HG11", "HG12", "HG13" ],
90 [ "HG1*", "HG11", "HG12", "HG13" ],
91 [ "HG2#", "HG21", "HG22", "HG23" ],
92 [ "HG2*", "HG21", "HG22", "HG23" ],
93 [ "HD#", "HD1", "HD2", "HD11", "HD12", "HD13", "HD21", "HD22", "HD23" ],
94 [ "HD*", "HD1", "HD2", "HD11", "HD12", "HD13", "HD21", "HD22", "HD23" ],
95 [ "HD1#", "HD11", "HD12" ],
96 [ "HD1*", "HD11", "HD12" ],
97 [ "HD2#", "HD21", "HD22" ],
98 [ "HD2*", "HD21", "HD22" ],
99 [ "HE#", "HE", "HE1", "HE2" ],
100 [ "HE*", "HE", "HE1", "HE2" ],
101 [ "HH#", "HH11", "HH12", "HH21", "HH22" ],
102 [ "HH*", "HH11", "HH12", "HH21", "HH22" ],
103 [ "HZ#", "HZ", "HZ1", "HZ2", "HZ3" ],
104 [ "HZ*", "HZ", "HZ1", "HZ2", "HZ3" ],
105 [ "HN", "H" ]
108 $ntranslated = 0;
109 $nuntransltd = 0;
110 sub transl_aname {
111 my $atom = shift(@_);
112 my $resnr = shift(@_);
114 if ( $trans == 1 ) {
115 for(my $i = 0; ($i < $ntbl); $i ++) {
116 if ($tblres[$i] eq $rname[$resnr]) {
117 if ($tblxplor[$i] eq $atom) {
118 $ntranslated++;
119 return $tblgmx[$i];
124 $nuntransltd++;
125 if ($debug > 1) {
126 printf "; No name change for $rname[$resnr] $atom\n";
128 return $atom;
131 sub expand_template {
132 my $atom = shift(@_);
133 my $rnum = shift(@_)-1;
134 my $bdone = 0;
135 my @atoms;
136 my $jj = 0;
138 for (my $tt=0; (($tt <= $#templates) && ($bdone == 0)); $tt++) {
139 $templ = $templates[$tt];
140 if ($atom eq $templ->[0]) {
141 for ($jj = 1; ($jj <= $#{$templ}); $jj++) {
142 push @atoms, transl_aname($templ->[$jj],$rnum);
144 $bdone = 1;
147 if ($bdone == 0) {
148 push @atoms, transl_aname($atom,$rnum);
150 if ($debug > 0) {
151 my $natom = $#atoms+1;
152 printf("; Found $natom elements for atom $rname[$rnum] %d $atom:",
153 $rnum+1);
154 for $aa ( @atoms ) {
155 printf " $aa";
157 printf "\n";
159 return @atoms;
162 if ($debug > 1) {
163 printf "; There are $#templates template entries\n";
164 for ($tt=0; ($tt <= $#templates); $tt++) {
165 $templ = $templates[$tt];
166 $ntempl = $#{$templ};
167 printf "; Item $tt ($templates[$tt][0]) has $ntempl entries\n";
171 # This index file holds numbers of restraints involving core atoms
172 @protcore = ( "H", "HN", "HA", "HA1", "HA2", "HB", "HB1", "HB2", "HB3", "HG", "HG1", "HG2", "HG3", "N", "O" );
173 open(CORE,">$core") || die "Can not open $core\n";
174 print CORE "[ core-restraints ]\n";
175 $ncore = 0;
177 $myindex = 0;
178 $linec = 0;
179 $npair = 0;
180 $nunresolved = 0;
181 while ($line = <STDIN>) {
182 @ttt = split('!',$line);
183 if ((index($ttt[0],"assign") >= 0) && (index($ttt[0],"!assign") < 0)) {
184 @tmp = split('\(',$ttt[0]);
185 # Find first argument
186 if (($rhaak = index($tmp[1],')')) < 0) {
187 printf "No ) in '$tmp[1]'\n";
189 $atres1 = substr($tmp[1],0,$rhaak);
190 @at1 = split('or',$atres1);
192 # Find second argument
193 if (($rhaak = index($tmp[2],')')) < 0) {
194 printf "No ) in '$tmp[2]'\n";
196 $atres2 = substr($tmp[2],0,$rhaak);
197 @at2 = split('or',$atres2);
199 @expdata = split('\)',$tmp[2]);
200 @dist = split(' ',$expdata[1]);
202 $bOK = 0;
203 $bCore = 1;
204 foreach $a1 ( @at1 ) {
205 @info1 = split(' ',$a1);
206 $r1 = $info1[1] - $res0;
207 @atoms1 = expand_template($info1[4],$r1);
209 foreach $a2 ( @at2 ) {
210 @info2 = split(' ',$a2);
211 $r2 = $info2[1] - $res0;
212 @atoms2 = expand_template($info2[4],$r2);
214 for ($i = 0; ($i < $natom) && ($resnr[$i] < $r1); $i++) { ; }
215 for ( ; ($i < $natom) && ($resnr[$i] == $r1); $i++) {
216 foreach $ii ( @atoms1 ) {
217 if ($ii eq $aname[$i]) {
218 $bCoreI = 0;
219 for $pp ( @protcore ) {
220 if ($ii eq $pp) {
221 $bCoreI = 1;
225 for ($j = 0; ($j < $natom) && ($resnr[$j] < $r2); $j++) { ; }
226 for ( ; ($j < $natom) && ($resnr[$j] == $r2); $j++) {
227 foreach $jj ( @atoms2 ) {
228 if ($jj eq $aname[$j]) {
229 $dd = 0.1*$dist[0];
230 $dminus = 0.1*$dist[1];
231 $dplus = 0.1*$dist[2];
232 $low = $dd-$dminus;
233 $up1 = $dd+$dplus;
234 $up2 = $up1+1;
235 printf("%5d %5d %5d %5d %5d %7.3f %7.3f %7.3f 1.0; res $r1 $ii - res $r2 $jj\n",$i+1,$j+1,1,$myindex,1,$low,$up1,$up2);
236 # Do some checks
237 $bOK = 1;
238 $bCoreJ = 0;
239 for $pp ( @protcore ) {
240 if ($jj eq $pp) {
241 $bCoreJ = 1;
244 if (($bCoreI == 0) || ($bCoreJ == 0)) {
245 if ($debug > 0) {
246 printf "; No core $ii ($bCoreI) $jj ($bCoreJ)\n";
248 $bCore = 0;
250 $npair++;
259 if ($bCore == 1) {
260 printf CORE "$myindex\n";
261 $ncore++;
263 if ($bOK == 0) {
264 printf "; UNRESOLVED: $ttt[0]\n";
265 $nunresolved++;
267 else {
268 $myindex++;
271 $linec++;
273 close CORE;
275 printf "; A total of $myindex lines with distance restraints were read from $linec input lines\n";
276 printf "; From this, $npair actual restraints were generated.\n";
277 printf "; A total of $nunresolved lines could not be interpreted\n";
278 printf "; In the process $ntranslated atoms had a name change\n";
279 printf "; However for $nuntransltd no name change could be found\n";
280 printf "; usually these are either HN->H renamings or not-existing atoms\n";
281 printf "; generated by the template expansion (e.g. HG#)\n";
282 printf "; A total of $ncore restraints were classified as core restraints\n";