2 C Read the PDB file and convert the peptide geometry into virtual-chain
9 include 'COMMON.INTERACT'
10 include 'COMMON.IOUNITS'
12 include 'COMMON.NAMES'
13 include 'COMMON.CONTROL'
15 include 'COMMON.SETUP'
16 include 'COMMON.SBRIDGE'
17 character*3 seq,atom,res
19 double precision sccor(3,50)
20 double precision e1(3),e2(3),e3(3)
21 integer rescode,iterter(maxres),cou
23 integer i,j,iii,ires,ires_old,ishift,ibeg
34 read (ipdbin,'(a80)',end=10) card
35 if (card(:5).eq.'HELIX') then
38 read(card(22:25),*) hfrag(1,nhfrag)
39 read(card(34:37),*) hfrag(2,nhfrag)
41 if (card(:5).eq.'SHEET') then
44 read(card(24:26),*) bfrag(1,nbfrag)
45 read(card(35:37),*) bfrag(2,nbfrag)
46 crc----------------------------------------
47 crc to be corrected !!!
48 bfrag(3,nbfrag)=bfrag(1,nbfrag)
49 bfrag(4,nbfrag)=bfrag(2,nbfrag)
50 crc----------------------------------------
52 if (card(:3).eq.'END') then
54 else if (card(:3).eq.'TER') then
57 itype(ires_old-1)=ntyp1
62 write (iout,*) "Chain ended",ires,ishift,ires_old
65 dc(j,ires)=sccor(j,iii)
68 call sccenter(ires,iii,sccor)
71 C Fish out the ATOM cards.
72 if (index(card(1:4),'ATOM').gt.0) then
73 read (card(14:16),'(a3)') atom
74 if (atom.eq.'CA' .or. atom.eq.'CH3') then
75 C Calculate the CM of the preceding residue.
79 dc(j,ires+nres)=sccor(j,iii)
82 call sccenter(ires,iii,sccor)
86 c write (iout,'(a80)') card
87 read (card(23:26),*) ires
88 read (card(18:20),'(a3)') res
91 if (res.ne.'GLY' .and. res.ne. 'ACE') then
95 c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
97 else if (ibeg.eq.2) then
99 ishift=-ires_old+ires-1
100 c write (iout,*) "New chain started",ires,ishift
104 c write (2,*) "ires",ires," ishift",ishift
105 if (res.eq.'ACE') then
108 itype(ires)=rescode(ires,res,0)
110 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
111 read(card(61:66),*) bfac(ires)
112 c if(me.eq.king.or..not.out1file)
113 c & write (iout,'(2i3,2x,a,3f8.3)')
114 c & ires,itype(ires),res,(c(j,ires),j=1,3)
117 sccor(j,iii)=c(j,ires)
119 else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and.
120 & atom.ne.'N ' .and. atom.ne.'C ') then
122 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
126 10 if(me.eq.king.or..not.out1file)
127 & write (iout,'(a,i5)') ' Nres: ',ires
128 C Calculate dummy residue coordinates inside the "chain" of a multichain
132 c write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
133 if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
134 if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
135 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
136 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
137 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
139 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
140 print *,i,'tu dochodze'
141 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
149 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
153 dcj=(c(j,i-2)-c(j,i-3))/2.0
154 if (dcj.eq.0) dcj=1.23591524223
159 else !itype(i+1).eq.ntyp1
161 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
162 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
169 c(j,i)=c(j,i+1)-1.9d0*e2(j)
173 dcj=(c(j,i+3)-c(j,i+2))/2.0
174 if (dcj.eq.0) dcj=1.23591524223
179 endif !itype(i+1).eq.ntyp1
180 endif !itype.eq.ntyp1
182 write (iout,*) "After loop in readpbd"
183 C Calculate the CM of the last side chain.
186 dc(j,ires)=sccor(j,iii)
189 call sccenter(ires,iii,sccor)
193 if (itype(nres).ne.10) then
197 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
198 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
205 c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
209 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
210 if (dcj.eq.0) dcj=1.23591524223
211 c(j,nres)=c(j,nres-1)+dcj
212 c(j,2*nres)=c(j,nres)
223 c(j,2*nres)=c(j,nres)
225 if (itype(1).eq.ntyp1) then
229 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
230 call refsys(2,3,4,e1,e2,e3,fail)
237 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
241 dcj=(c(j,4)-c(j,3))/2.0
247 C Calculate internal coordinates.
248 if(me.eq.king.or..not.out1file)then
250 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
251 & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
252 & (c(j,nres+ires),j=1,3)
256 c write(iout,*)"before int_from_cart nres",nres
257 call int_from_cart(.true.,.false.)
264 dc(j,i)=c(j,i+1)-c(j,i)
265 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
267 c write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
272 dc(j,i+nres)=c(j,i+nres)-c(j,i)
273 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
275 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
278 call sc_loc_geom(.false.)
279 call int_from_cart1(.false.)
281 C Copy the coordinates to reference coordinates
285 cref(j,i+nres)=c(j,i+nres)
288 100 format (//' alpha-carbon coordinates ',
289 & ' centroid coordinates'/
290 1 ' ', 6X,'X',11X,'Y',11X,'Z',
291 & 10X,'X',11X,'Y',11X,'Z')
292 110 format (a,'(',i3,')',6f12.5)
296 bfrag(i,j)=bfrag(i,j)-ishift
302 hfrag(i,j)=hfrag(i,j)-ishift
307 c---------------------------------------------------------------------------
308 subroutine readpdb_template(k)
309 C Read the PDB file for read_constr_homology with read2sigma
310 C and convert the peptide geometry into virtual-chain geometry.
313 include 'COMMON.LOCAL'
315 include 'COMMON.CHAIN'
316 include 'COMMON.INTERACT'
317 include 'COMMON.IOUNITS'
319 include 'COMMON.NAMES'
320 include 'COMMON.CONTROL'
321 include 'COMMON.FRAG'
322 include 'COMMON.SETUP'
323 integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
325 logical lprn /.false./,fail
326 double precision e1(3),e2(3),e3(3)
327 double precision dcj,efree_temp
331 double precision sccor(3,20)
332 integer rescode,iterter(maxres)
339 c write (2,*) "UNRES_PDB",unres_pdb
347 read (ipdbin,'(a80)',end=10) card
348 if (card(:3).eq.'END') then
350 else if (card(:3).eq.'TER') then
353 itype(ires_old-1)=ntyp1
354 iterter(ires_old-1)=1
355 itype(ires_old)=ntyp1
358 c write (iout,*) "Chain ended",ires,ishift,ires_old
361 dc(j,ires)=sccor(j,iii)
364 call sccenter(ires,iii,sccor)
367 C Fish out the ATOM cards.
368 if (index(card(1:4),'ATOM').gt.0) then
369 read (card(12:16),*) atom
370 c write (iout,*) "! ",atom," !",ires
371 c if (atom.eq.'CA' .or. atom.eq.'CH3') then
372 read (card(23:26),*) ires
373 read (card(18:20),'(a3)') res
374 c write (iout,*) "ires",ires,ires-ishift+ishift1,
375 c & " ires_old",ires_old
376 c write (iout,*) "ishift",ishift," ishift1",ishift1
377 c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
378 if (ires-ishift+ishift1.ne.ires_old) then
379 C Calculate the CM of the preceding residue.
383 dc(j,ires)=sccor(j,iii)
386 call sccenter(ires_old,iii,sccor)
391 if (res.eq.'Cl-' .or. res.eq.'Na+') then
394 else if (ibeg.eq.1) then
395 c write (iout,*) "BEG ires",ires
397 if (res.ne.'GLY' .and. res.ne. 'ACE') then
401 ires=ires-ishift+ishift1
403 c write (iout,*) "ishift",ishift," ires",ires,
404 c & " ires_old",ires_old
405 c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
407 else if (ibeg.eq.2) then
409 ishift=-ires_old+ires-1
411 c write (iout,*) "New chain started",ires,ishift
414 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
415 ires=ires-ishift+ishift1
418 if (res.eq.'ACE' .or. res.eq.'NHE') then
421 itype(ires)=rescode(ires,res,0)
424 ires=ires-ishift+ishift1
426 c write (iout,*) "ires_old",ires_old," ires",ires
427 c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
430 c write (2,*) "ires",ires," res ",res," ity",ity
431 if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
432 & res.eq.'NHE'.and.atom(:2).eq.'HN') then
433 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
434 c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
436 write (iout,'(2i3,2x,a,3f8.3)')
437 & ires,itype(ires),res,(c(j,ires),j=1,3)
441 sccor(j,iii)=c(j,ires)
443 if (ishift.ne.0) then
444 ires_ca=ires+ishift-ishift1
448 c write (*,*) card(23:27),ires,itype(ires)
449 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
450 & atom.ne.'N' .and. atom.ne.'C' .and.
451 & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
452 & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
453 c write (iout,*) "sidechain ",atom
455 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
459 10 if(me.eq.king.or..not.out1file)
460 & write (iout,'(a,i5)') ' Nres: ',ires
461 C Calculate dummy residue coordinates inside the "chain" of a multichain
465 c write (iout,*) i,itype(i),itype(i+1)
466 if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
467 if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
468 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
469 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
470 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
472 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
473 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
480 c(j,i)=c(j,i-1)-1.9d0*e2(j)
484 dcj=(c(j,i-2)-c(j,i-3))/2.0
485 if (dcj.eq.0) dcj=1.23591524223
490 else !itype(i+1).eq.ntyp1
492 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
493 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
500 c(j,i)=c(j,i+1)-1.9d0*e2(j)
504 dcj=(c(j,i+3)-c(j,i+2))/2.0
505 if (dcj.eq.0) dcj=1.23591524223
510 endif !itype(i+1).eq.ntyp1
511 endif !itype.eq.ntyp1
513 C Calculate the CM of the last side chain.
516 dc(j,ires)=sccor(j,iii)
519 call sccenter(ires,iii,sccor)
523 if (itype(nres).ne.10) then
527 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
528 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
535 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
539 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
540 if (dcj.eq.0) dcj=1.23591524223
541 c(j,nres)=c(j,nres-1)+dcj
542 c(j,2*nres)=c(j,nres)
553 c(j,2*nres)=c(j,nres)
555 if (itype(1).eq.ntyp1) then
559 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
560 call refsys(2,3,4,e1,e2,e3,fail)
567 c(j,1)=c(j,2)-1.9d0*e2(j)
571 dcj=(c(j,4)-c(j,3))/2.0
577 C Copy the coordinates to reference coordinates
583 C Calculate internal coordinates.
584 if (out_template_coord) then
586 & "Cartesian coordinates of the reference structure"
587 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
588 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
590 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
591 & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
592 & (c(j,ires+nres),j=1,3)
595 C Calculate internal coordinates.
596 call int_from_cart(.true.,.true.)
597 call sc_loc_geom(.false.)
604 dc(j,i)=c(j,i+1)-c(j,i)
605 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
610 dc(j,i+nres)=c(j,i+nres)-c(j,i)
611 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
613 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
619 cref(j,i+nres)=c(j,i+nres)