2 C Read the PDB file and convert the peptide geometry into virtual-chain
4 implicit real*8 (a-h,o-z)
9 include 'COMMON.INTERACT'
10 include 'COMMON.IOUNITS'
12 include 'COMMON.NAMES'
13 include 'COMMON.CONTROL'
14 include 'COMMON.DISTFIT'
15 include 'COMMON.SETUP'
17 integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
19 logical lprn /.false./,fail
20 double precision e1(3),e2(3),e3(3)
21 double precision dcj,efree_temp
25 double precision sccor(3,50)
31 c write (2,*) "UNRES_PDB",unres_pdb
39 read (ipdbin,'(a80)',end=10) card
40 c write (iout,'(a)') card
41 if (card(:5).eq.'HELIX') then
44 read(card(22:25),*) hfrag(1,nhfrag)
45 read(card(34:37),*) hfrag(2,nhfrag)
47 if (card(:5).eq.'SHEET') then
50 read(card(24:26),*) bfrag(1,nbfrag)
51 read(card(35:37),*) bfrag(2,nbfrag)
52 crc----------------------------------------
53 crc to be corrected !!!
54 bfrag(3,nbfrag)=bfrag(1,nbfrag)
55 bfrag(4,nbfrag)=bfrag(2,nbfrag)
56 crc----------------------------------------
58 if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
60 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
61 C Fish out the ATOM cards.
62 if (index(card(1:4),'ATOM').gt.0) then
63 read (card(12:16),*) atom
64 c write (iout,*) "! ",atom," !",ires
65 c if (atom.eq.'CA' .or. atom.eq.'CH3') then
66 read (card(23:26),*) ires
67 read (card(18:20),'(a3)') res
68 c write (iout,*) "ires",ires,ires-ishift+ishift1,
69 c & " ires_old",ires_old
70 c write (iout,*) "ishift",ishift," ishift1",ishift1
71 c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
72 if (ires-ishift+ishift1.ne.ires_old) then
73 C Calculate the CM of the preceding residue.
74 c if (ibeg.eq.0) call sccenter(ires,iii,sccor)
76 c write (iout,*) "Calculating sidechain center iii",iii
79 dc(j,ires)=sccor(j,iii)
82 call sccenter(ires_old,iii,sccor)
87 if (res.eq.'Cl-' .or. res.eq.'Na+') then
90 else if (ibeg.eq.1) then
91 c write (iout,*) "BEG ires",ires
93 if (res.ne.'GLY' .and. res.ne. 'ACE') then
97 ires=ires-ishift+ishift1
99 c write (iout,*) "ishift",ishift," ires",ires,
100 c & " ires_old",ires_old
103 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
104 ires=ires-ishift+ishift1
107 if (res.eq.'ACE' .or. res.eq.'NHE') then
110 itype(ires)=rescode(ires,res,0)
113 ires=ires-ishift+ishift1
115 c write (iout,*) "ires_old",ires_old," ires",ires
116 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
119 c write (2,*) "ires",ires," res ",res," ity",ity
120 if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
121 & res.eq.'NHE'.and.atom(:2).eq.'HN') then
122 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
123 c write (iout,*) "backbone ",atom
125 write (iout,'(2i3,2x,a,3f8.3)')
126 & ires,itype(ires),res,(c(j,ires),j=1,3)
130 sccor(j,iii)=c(j,ires)
132 if (ishift.ne.0) then
133 ires_ca=ires+ishift-ishift1
137 c write (*,*) card(23:27),ires,itype(ires)
138 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
139 & atom.ne.'N' .and. atom.ne.'C' .and.
140 & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
141 & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
142 c write (iout,*) "sidechain ",atom
144 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
150 write (iout,'(a,i5)') ' Number of residues found: ',ires
152 if (ires.eq.0) return
153 C Calculate the CM of the last side chain.
157 dc(j,ires)=sccor(j,iii)
160 call sccenter(ires,iii,sccor)
166 if (itype(nres).ne.10) then
170 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
171 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
178 c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
182 dcj=c(j,nres-2)-c(j,nres-3)
183 c(j,nres)=c(j,nres-1)+dcj
184 c(j,2*nres)=c(j,nres)
195 c(j,2*nres)=c(j,nres)
197 if (itype(1).eq.ntyp1) then
201 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
202 call refsys(2,3,4,e1,e2,e3,fail)
209 c(j,1)=c(j,2)-3.8d0*e2(j)
219 C Copy the coordinates to reference coordinates
225 C Calculate internal coordinates.
228 & "Cartesian coordinates of the reference structure"
229 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
230 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
232 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
233 & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
234 & (c(j,ires+nres),j=1,3)
237 C Calculate internal coordinates.
238 if(me.eq.king.or..not.out1file)then
240 & "Backbone and SC coordinates as read from the PDB"
242 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
243 & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
244 & (c(j,nres+ires),j=1,3)
247 call int_from_cart(.true.,.false.)
248 call sc_loc_geom(.false.)
255 dc(j,i)=c(j,i+1)-c(j,i)
256 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
261 dc(j,i+nres)=c(j,i+nres)-c(j,i)
262 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
264 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
268 C Copy the coordinates to reference coordinates
278 bfrag(i,j)=bfrag(i,j)-ishift
284 hfrag(i,j)=hfrag(i,j)-ishift
290 c-----------------------------------------------------------------------
291 subroutine readpdb_template(k)
292 C Read the PDB file with gaps for read_constr_homology with read2sigma
293 C and convert the peptide geometry into virtual-chain geometry.
294 implicit real*8 (a-h,o-z)
296 include 'COMMON.LOCAL'
298 include 'COMMON.CHAIN'
299 include 'COMMON.INTERACT'
300 include 'COMMON.IOUNITS'
302 include 'COMMON.NAMES'
303 include 'COMMON.CONTROL'
304 include 'COMMON.DISTFIT'
305 include 'COMMON.SETUP'
307 integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
309 logical lprn /.false./,fail
310 double precision e1(3),e2(3),e3(3)
311 double precision dcj,efree_temp
315 double precision sccor(3,50)
321 c write (2,*) "UNRES_PDB",unres_pdb
329 read (ipdbin,'(a80)',end=10) card
330 c write (iout,'(a)') card
331 if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
332 C Fish out the ATOM cards.
333 if (index(card(1:4),'ATOM').gt.0) then
334 read (card(12:16),*) atom
335 c write (iout,*) "! ",atom," !",ires
336 c if (atom.eq.'CA' .or. atom.eq.'CH3') then
337 read (card(23:26),*) ires
338 read (card(18:20),'(a3)') res
339 c write (iout,*) "ires",ires,ires-ishift+ishift1,
340 c & " ires_old",ires_old
341 c write (iout,*) "ishift",ishift," ishift1",ishift1
342 c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
343 if (ires-ishift+ishift1.ne.ires_old) then
344 C Calculate the CM of the preceding residue.
345 c if (ibeg.eq.0) call sccenter(ires,iii,sccor)
347 c write (iout,*) "Calculating sidechain center iii",iii
350 dc(j,ires)=sccor(j,iii)
353 call sccenter(ires_old,iii,sccor)
358 if (res.eq.'Cl-' .or. res.eq.'Na+') then
361 else if (ibeg.eq.1) then
362 c write (iout,*) "BEG ires",ires
364 if (res.ne.'GLY' .and. res.ne. 'ACE') then
368 ires=ires-ishift+ishift1
370 c write (iout,*) "ishift",ishift," ires",ires,
371 c & " ires_old",ires_old
374 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
375 ires=ires-ishift+ishift1
378 if (res.eq.'ACE' .or. res.eq.'NHE') then
381 itype(ires)=rescode(ires,res,0)
384 ires=ires-ishift+ishift1
386 c write (iout,*) "ires_old",ires_old," ires",ires
387 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
390 c write (2,*) "ires",ires," res ",res," ity",ity
391 if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
392 & res.eq.'NHE'.and.atom(:2).eq.'HN') then
393 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
394 c write (iout,*) "backbone ",atom
396 write (iout,'(2i3,2x,a,3f8.3)')
397 & ires,itype(ires),res,(c(j,ires),j=1,3)
401 sccor(j,iii)=c(j,ires)
403 if (ishift.ne.0) then
404 ires_ca=ires+ishift-ishift1
408 c write (*,*) card(23:27),ires,itype(ires)
409 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
410 & atom.ne.'N' .and. atom.ne.'C' .and.
411 & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
412 & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
413 c write (iout,*) "sidechain ",atom
415 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
421 write (iout,'(a,i5)') ' Number of residues found: ',ires
423 if (ires.eq.0) return
424 C Calculate the CM of the last side chain.
428 dc(j,ires)=sccor(j,iii)
431 call sccenter(ires,iii,sccor)
437 if (itype(nres).ne.10) then
441 dcj=c(j,nres-2)-c(j,nres-3)
442 c(j,nres)=c(j,nres-1)+dcj
443 c(j,2*nres)=c(j,nres)
453 c(j,2*nres)=c(j,nres)
455 if (itype(1).eq.ntyp1) then
459 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
460 call refsys(2,3,4,e1,e2,e3,fail)
467 c(j,1)=c(j,2)-3.8d0*e2(j)
477 C Calculate internal coordinates.
480 & "Cartesian coordinates of the reference structure"
481 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
482 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
484 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
485 & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
486 & (c(j,ires+nres),j=1,3)
489 C Calculate internal coordinates.
490 if(me.eq.king.or..not.out1file)then
492 & "Backbone and SC coordinates as read from the PDB"
494 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
495 & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
496 & (c(j,nres+ires),j=1,3)
499 call int_from_cart(.true.,.false.)
500 call sc_loc_geom(.false.)
507 dc(j,i)=c(j,i+1)-c(j,i)
508 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
513 dc(j,i+nres)=c(j,i+nres)-c(j,i)
514 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
516 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
520 C Copy the coordinates to reference coordinates