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,ishift1,ibeg
24 double precision dcj,efree_temp
35 read (ipdbin,'(a80)',end=10) card
36 ! write (iout,'(a)') card
37 if (card(:5).eq.'HELIX') then
40 read(card(22:25),*) hfrag(1,nhfrag)
41 read(card(34:37),*) hfrag(2,nhfrag)
43 if (card(:5).eq.'SHEET') then
46 read(card(24:26),*) bfrag(1,nbfrag)
47 read(card(35:37),*) bfrag(2,nbfrag)
48 !rc----------------------------------------
49 !rc to be corrected !!!
50 bfrag(3,nbfrag)=bfrag(1,nbfrag)
51 bfrag(4,nbfrag)=bfrag(2,nbfrag)
52 !rc----------------------------------------
54 if (card(:3).eq.'END') then
56 else if (card(:3).eq.'TER') then
59 itype(ires_old-1)=ntyp1
64 ! write (iout,*) "Chain ended",ires,ishift,ires_old
67 dc(j,ires)=sccor(j,iii)
70 call sccenter(ires,iii,sccor)
75 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
76 ! Fish out the ATOM cards.
77 if (index(card(1:4),'ATOM').gt.0) then
78 read (card(12:16),*) atom
79 c write (2,'(a)') card
80 ! write (iout,*) "! ",atom," !",ires
81 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
82 read (card(23:26),*) ires
83 read (card(18:20),'(a3)') res
84 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
85 ! & " ires_old",ires_old
86 ! write (iout,*) "ishift",ishift," ishift1",ishift1
87 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
88 if (ires-ishift+ishift1.ne.ires_old) then
89 ! Calculate the CM of the preceding residue.
90 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
92 ! write (iout,*) "Calculating sidechain center iii",iii
95 dc(j,ires+nres)=sccor(j,iii)
98 call sccenter(ires_old,iii,sccor)
103 if (res.eq.'Cl-' .or. res.eq.'Na+') then
106 else if (ibeg.eq.1) then
107 c write (iout,*) "BEG ires",ires
109 if (res.ne.'GLY' .and. res.ne. 'ACE') then
113 ires=ires-ishift+ishift1
115 ! write (iout,*) "ishift",ishift," ires",ires,&
116 ! " ires_old",ires_old
118 else if (ibeg.eq.2) then
120 ishift=-ires_old+ires-1 !!!!!
121 ishift1=ishift1-1 !!!!!
122 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
123 ires=ires-ishift+ishift1
127 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
128 ires=ires-ishift+ishift1
131 if (res.eq.'ACE' .or. res.eq.'NHE') then
134 itype(ires)=rescode(ires,res,0)
137 ires=ires-ishift+ishift1
139 ! write (iout,*) "ires_old",ires_old," ires",ires
140 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
143 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
144 if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
145 & res.eq.'NHE'.and.atom(:2).eq.'HN') then
146 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
147 ! write (iout,*) "backbone ",atom
149 write (iout,'(2i3,2x,a,3f8.3)')
150 & ires,itype(ires),res,(c(j,ires),j=1,3)
154 sccor(j,iii)=c(j,ires)
156 c write (2,*) card(23:27),ires,itype(ires),iii
157 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
158 & atom.ne.'N' .and. atom.ne.'C' .and.
159 & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
160 & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
161 ! write (iout,*) "sidechain ",atom
163 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
164 c write (2,*) "iii",iii
168 10 if(me.eq.king.or..not.out1file)
169 & write (iout,'(a,i5)') ' Nres: ',ires
170 c write (iout,*) "iii",iii
171 C Calculate dummy residue coordinates inside the "chain" of a multichain
175 c write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
176 if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
177 if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
178 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
179 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
180 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
182 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
183 print *,i,'tu dochodze'
184 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
192 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
196 dcj=(c(j,i-2)-c(j,i-3))/2.0
197 if (dcj.eq.0) dcj=1.23591524223
202 else !itype(i+1).eq.ntyp1
204 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
205 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
212 c(j,i)=c(j,i+1)-1.9d0*e2(j)
216 dcj=(c(j,i+3)-c(j,i+2))/2.0
217 if (dcj.eq.0) dcj=1.23591524223
222 endif !itype(i+1).eq.ntyp1
223 endif !itype.eq.ntyp1
225 write (iout,*) "After loop in readpbd"
226 C Calculate the CM of the last side chain.
229 dc(j,ires)=sccor(j,iii)
232 c write (iout,*) "Calling sccenter iii",iii
233 call sccenter(ires,iii,sccor)
237 if (itype(nres).ne.10) then
241 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
242 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
249 c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
253 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
254 if (dcj.eq.0) dcj=1.23591524223
255 c(j,nres)=c(j,nres-1)+dcj
256 c(j,2*nres)=c(j,nres)
267 c(j,2*nres)=c(j,nres)
269 if (itype(1).eq.ntyp1) then
273 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
274 call refsys(2,3,4,e1,e2,e3,fail)
281 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
285 dcj=(c(j,4)-c(j,3))/2.0
291 C Calculate internal coordinates.
292 if(me.eq.king.or..not.out1file)then
294 & "Cartesian coordinates of the reference structure"
295 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
296 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
298 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
299 & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
300 & (c(j,ires+nres),j=1,3)
304 c write(iout,*)"before int_from_cart nres",nres
305 call int_from_cart(.true.,.false.)
312 dc(j,i)=c(j,i+1)-c(j,i)
313 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
315 c write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
320 dc(j,i+nres)=c(j,i+nres)-c(j,i)
321 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
323 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
326 call sc_loc_geom(.false.)
327 call int_from_cart1(.false.)
329 C Copy the coordinates to reference coordinates
333 cref(j,i+nres)=c(j,i+nres)
336 100 format (//' alpha-carbon coordinates ',
337 & ' centroid coordinates'/
338 1 ' ', 6X,'X',11X,'Y',11X,'Z',
339 & 10X,'X',11X,'Y',11X,'Z')
340 110 format (a,'(',i3,')',6f12.5)
344 bfrag(i,j)=bfrag(i,j)-ishift
350 hfrag(i,j)=hfrag(i,j)-ishift
355 c---------------------------------------------------------------------------
356 subroutine readpdb_template(k)
357 C Read the PDB file for read_constr_homology with read2sigma
358 C and convert the peptide geometry into virtual-chain geometry.
361 include 'COMMON.LOCAL'
363 include 'COMMON.CHAIN'
364 include 'COMMON.INTERACT'
365 include 'COMMON.IOUNITS'
367 include 'COMMON.NAMES'
368 include 'COMMON.CONTROL'
369 include 'COMMON.FRAG'
370 include 'COMMON.SETUP'
371 integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
373 logical lprn /.false./,fail
374 double precision e1(3),e2(3),e3(3)
375 double precision dcj,efree_temp
379 double precision sccor(3,20)
380 integer rescode,iterter(maxres)
387 c write (2,*) "UNRES_PDB",unres_pdb
395 read (ipdbin,'(a80)',end=10) card
396 if (card(:3).eq.'END') then
398 else if (card(:3).eq.'TER') then
401 itype(ires_old-1)=ntyp1
402 iterter(ires_old-1)=1
403 itype(ires_old)=ntyp1
406 c write (iout,*) "Chain ended",ires,ishift,ires_old
409 dc(j,ires)=sccor(j,iii)
412 call sccenter(ires,iii,sccor)
415 C Fish out the ATOM cards.
416 if (index(card(1:4),'ATOM').gt.0) then
417 read (card(12:16),*) atom
418 c write (iout,*) "! ",atom," !",ires
419 c if (atom.eq.'CA' .or. atom.eq.'CH3') then
420 read (card(23:26),*) ires
421 read (card(18:20),'(a3)') res
422 c write (iout,*) "ires",ires,ires-ishift+ishift1,
423 c & " ires_old",ires_old
424 c write (iout,*) "ishift",ishift," ishift1",ishift1
425 c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
426 if (ires-ishift+ishift1.ne.ires_old) then
427 C Calculate the CM of the preceding residue.
431 dc(j,ires)=sccor(j,iii)
434 call sccenter(ires_old,iii,sccor)
439 if (res.eq.'Cl-' .or. res.eq.'Na+') then
442 else if (ibeg.eq.1) then
443 c write (iout,*) "BEG ires",ires
445 if (res.ne.'GLY' .and. res.ne. 'ACE') then
449 ires=ires-ishift+ishift1
451 c write (iout,*) "ishift",ishift," ires",ires,
452 c & " ires_old",ires_old
453 c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
455 else if (ibeg.eq.2) then
457 ishift=-ires_old+ires-1
459 c write (iout,*) "New chain started",ires,ishift
462 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
463 ires=ires-ishift+ishift1
466 if (res.eq.'ACE' .or. res.eq.'NHE') then
469 itype(ires)=rescode(ires,res,0)
472 ires=ires-ishift+ishift1
474 c write (iout,*) "ires_old",ires_old," ires",ires
475 c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
478 c write (2,*) "ires",ires," res ",res," ity",ity
479 if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
480 & res.eq.'NHE'.and.atom(:2).eq.'HN') then
481 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
482 c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
484 write (iout,'(2i3,2x,a,3f8.3)')
485 & ires,itype(ires),res,(c(j,ires),j=1,3)
489 sccor(j,iii)=c(j,ires)
491 if (ishift.ne.0) then
492 ires_ca=ires+ishift-ishift1
496 c write (*,*) card(23:27),ires,itype(ires)
497 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
498 & atom.ne.'N' .and. atom.ne.'C' .and.
499 & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
500 & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
501 c write (iout,*) "sidechain ",atom
503 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
507 10 if(me.eq.king.or..not.out1file)
508 & write (iout,'(a,i5)') ' Nres: ',ires
509 C Calculate dummy residue coordinates inside the "chain" of a multichain
513 c write (iout,*) i,itype(i),itype(i+1)
514 if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
515 if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
516 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
517 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
518 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
520 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
521 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
528 c(j,i)=c(j,i-1)-1.9d0*e2(j)
532 dcj=(c(j,i-2)-c(j,i-3))/2.0
533 if (dcj.eq.0) dcj=1.23591524223
538 else !itype(i+1).eq.ntyp1
540 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
541 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
548 c(j,i)=c(j,i+1)-1.9d0*e2(j)
552 dcj=(c(j,i+3)-c(j,i+2))/2.0
553 if (dcj.eq.0) dcj=1.23591524223
558 endif !itype(i+1).eq.ntyp1
559 endif !itype.eq.ntyp1
561 C Calculate the CM of the last side chain.
564 dc(j,ires)=sccor(j,iii)
567 call sccenter(ires,iii,sccor)
571 if (itype(nres).ne.10) then
575 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
576 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
583 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
587 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
588 if (dcj.eq.0) dcj=1.23591524223
589 c(j,nres)=c(j,nres-1)+dcj
590 c(j,2*nres)=c(j,nres)
601 c(j,2*nres)=c(j,nres)
603 if (itype(1).eq.ntyp1) then
607 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
608 call refsys(2,3,4,e1,e2,e3,fail)
615 c(j,1)=c(j,2)-1.9d0*e2(j)
619 dcj=(c(j,4)-c(j,3))/2.0
625 C Copy the coordinates to reference coordinates
631 C Calculate internal coordinates.
632 if (out_template_coord) then
634 & "Cartesian coordinates of the reference structure"
635 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
636 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
638 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
639 & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
640 & (c(j,ires+nres),j=1,3)
643 C Calculate internal coordinates.
644 call int_from_cart(.true.,.true.)
645 call sc_loc_geom(.false.)
652 dc(j,i)=c(j,i+1)-c(j,i)
653 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
658 dc(j,i+nres)=c(j,i+nres)-c(j,i)
659 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
661 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
667 cref(j,i+nres)=c(j,i+nres)