2 C Read the PDB file and convert the peptide geometry into virtual-chain
6 include 'DIMENSIONS.ZSCOPT'
7 include 'COMMON.CONTROL'
10 include 'COMMON.CHAIN'
11 include 'COMMON.INTERACT'
12 include 'COMMON.IOUNITS'
14 include 'COMMON.NAMES'
15 include 'COMMON.SBRIDGE'
17 character*3 seq,atom,res
19 double precision sccor(3,50)
20 integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
22 integer rescode,kkk,lll,icha,cou,kupa,iprzes
23 logical lsecondary,sccalc
24 integer iterter(maxres)
25 double precision efree_temp
32 read (ipdbin,'(a80)',end=10) card
33 c write (iout,'(a)') card
34 if (card(:5).eq.'HELIX') then
37 read(card(22:25),*) hfrag(1,nhfrag)
38 read(card(34:37),*) hfrag(2,nhfrag)
40 if (card(:5).eq.'SHEET') then
43 read(card(24:26),*) bfrag(1,nbfrag)
44 read(card(35:37),*) bfrag(2,nbfrag)
45 !rc----------------------------------------
46 !rc to be corrected !!!
47 bfrag(3,nbfrag)=bfrag(1,nbfrag)
48 bfrag(4,nbfrag)=bfrag(2,nbfrag)
49 !rc----------------------------------------
51 if (card(:3).eq.'END') then
53 else if (card(:3).eq.'TER') then
56 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)
74 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
75 ! Fish out the ATOM cards.
76 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)
104 if (res.eq.'Cl-' .or. res.eq.'Na+') then
107 else if (ibeg.eq.1) then
108 c write (iout,*) "BEG ires",ires
110 if (res.ne.'GLY' .and. res.ne. 'ACE') then
114 ires=ires-ishift+ishift1
116 ! write (iout,*) "ishift",ishift," ires",ires,&
117 ! " ires_old",ires_old
119 else if (ibeg.eq.2) then
121 ishift=-ires_old+ires-1 !!!!!
122 ishift1=ishift1-1 !!!!!
123 write (iout,*) "New chain started",ires,ishift,ishift1,"!"
124 ires=ires-ishift+ishift1
128 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
129 ires=ires-ishift+ishift1
132 if (res.eq.'ACE' .or. res.eq.'NHE') then
135 itype(ires)=rescode(ires,res,0)
138 ires=ires-ishift+ishift1
140 ! write (iout,*) "ires_old",ires_old," ires",ires
141 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
144 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
145 if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
146 & res.eq.'NHE'.and.atom(:2).eq.'HN') then
147 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
148 ! write (iout,*) "backbone ",atom
150 write (iout,'(2i3,2x,a,3f8.3)')
151 & ires,itype(ires),res,(c(j,ires),j=1,3)
155 sccor(j,iii)=c(j,ires)
157 c write (2,*) card(23:27),ires,itype(ires),iii
158 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
159 & atom.ne.'N' .and. atom.ne.'C' .and.
160 & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
161 & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
162 ! write (iout,*) "sidechain ",atom
164 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
165 c write (2,*) "iii",iii
169 10 write (iout,'(a,i5)') ' Nres: ',ires
170 C Calculate dummy residue coordinates inside the "chain" of a multichain
174 c write (iout,*) i,itype(i)
176 if (itype(i).eq.ntyp1) then
177 if (itype(i+1).eq.ntyp1) 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
181 C if (unres_pdb) then
182 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
183 C call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
190 C c(j,i)=c(j,i-1)-1.9d0*e2(j)
194 dcj=(c(j,i-2)-c(j,i-3))/2.0
199 else !itype(i+1).eq.ntyp1
200 C if (unres_pdb) then
201 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
202 C call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
209 C c(j,i)=c(j,i+1)-1.9d0*e2(j)
213 dcj=(c(j,i+3)-c(j,i+2))/2.0
218 endif !itype(i+1).eq.ntyp1
219 endif !itype.eq.ntyp1
221 C Calculate the CM of the last side chain.
222 if (.not. sccalc) call sccenter(ires,iii,sccor)
225 if (itype(nres).ne.10) then
229 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
230 c(j,nres)=c(j,nres-1)+dcj
231 c(j,2*nres)=c(j,nres)
241 c(j,2*nres)=c(j,nres)
243 if (itype(1).eq.ntyp1) then
247 dcj=(c(j,4)-c(j,3))/2.0
252 C Calculate internal coordinates.
255 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
256 & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
257 & (c(j,nres+ires),j=1,3)
259 call int_from_cart(.true.,.false.)
263 dc(j,i)=c(j,i+1)-c(j,i)
264 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
269 dc(j,i+nres)=c(j,i+nres)-c(j,i)
270 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
272 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
276 C Copy the coordinates to reference coordinates
280 cref(j,i+nres)=c(j,i+nres)
283 100 format ('Residue alpha-carbon coordinates ',
284 & ' centroid coordinates'/
285 1 ' ', 6X,'X',7X,'Y',7X,'Z',
286 & 12X,'X',7X,'Y',7X,'Z')
287 110 format (a,'(',i3,')',6f12.5)
292 c---------------------------------------------------------------------------
293 subroutine int_from_cart(lside,lprn)
296 include 'DIMENSIONS.ZSCOPT'
297 include 'COMMON.LOCAL'
299 include 'COMMON.CHAIN'
300 include 'COMMON.INTERACT'
301 include 'COMMON.IOUNITS'
303 include 'COMMON.NAMES'
304 character*3 seq,atom,res
306 double precision sccor(3,50)
308 double precision dist,alpha,beta,di
313 & 'Internal coordinates calculated from crystal structure.'
315 write (iout,'(8a)') ' Res ',' dvb',' Theta',
316 & ' Phi',' Dsc_id',' Dsc',' Alpha',
319 write (iout,'(4a)') ' Res ',' dvb',' Theta',
325 c write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
326 if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
327 & (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
328 write (iout,'(a,i4)') 'Bad Cartesians for residue',i
332 vbld_inv(i)=1.0d0/vbld(i)
333 theta(i+1)=alpha(i-1,i,i+1)
334 if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
336 c if (itype(1).eq.ntyp1) then
338 c c(j,1)=c(j,2)+(c(j,3)-c(j,4))
341 c if (itype(nres).eq.ntyp1) then
343 c c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
349 c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
354 if (itype(i).ne.10) then
355 vbld_inv(i+nres)=1.0d0/di
357 vbld_inv(i+nres)=0.0d0
360 alph(i)=alpha(nres+i,i,maxres2)
361 omeg(i)=beta(nres+i,i,maxres2,i+1)
364 alph(i)=alpha(nres+i,i,maxres2)
365 omeg(i)=beta(nres+i,i,maxres2,i+1)
368 & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
369 & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
370 & rad2deg*alph(i),rad2deg*omeg(i)
375 write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
376 & rad2deg*theta(i),rad2deg*phi(i)
381 c---------------------------------------------------------------------------
382 subroutine sccenter(ires,nscat,sccor)
385 include 'COMMON.CHAIN'
386 integer ires,nscat,i,j
387 double precision sccor(3,50),sccmj
391 sccmj=sccmj+sccor(j,i)
393 dc(j,ires)=sccmj/nscat
397 c---------------------------------------------------------------------------
398 subroutine sc_loc_geom(lprn)
399 implicit real*8 (a-h,o-z)
401 include 'DIMENSIONS.ZSCOPT'
402 include 'COMMON.LOCAL'
404 include 'COMMON.CHAIN'
405 include 'COMMON.INTERACT'
406 include 'COMMON.IOUNITS'
408 include 'COMMON.NAMES'
409 include 'COMMON.CONTROL'
410 include 'COMMON.SETUP'
411 double precision x_prime(3),y_prime(3),z_prime(3)
415 dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
419 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
421 dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
425 dc_norm(j,i+nres)=0.0d0
430 costtab(i+1) =dcos(theta(i+1))
431 sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
432 cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
433 sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
434 cosfac2=0.5d0/(1.0d0+costtab(i+1))
435 cosfac=dsqrt(cosfac2)
436 sinfac2=0.5d0/(1.0d0-costtab(i+1))
437 sinfac=dsqrt(sinfac2)
439 if (it.ne.10 .and. itype(i).ne.ntyp1) then
441 C Compute the axes of tghe local cartesian coordinates system; store in
442 c x_prime, y_prime and z_prime
450 x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
451 y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
453 c write (iout,*) "x_prime",(x_prime(j),j=1,3)
454 c write (iout,*) "y_prime",(y_prime(j),j=1,3)
455 call vecpr(x_prime,y_prime,z_prime)
456 c write (iout,*) "z_prime",(z_prime(j),j=1,3)
458 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
459 C to local coordinate system. Store in xx, yy, zz.
465 xx = xx + x_prime(j)*dc_norm(j,i+nres)
466 yy = yy + y_prime(j)*dc_norm(j,i+nres)
467 zz = zz + z_prime(j)*dc_norm(j,i+nres)
480 write (iout,*) "xxref,yyref,zzref"
483 write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
489 c---------------------------------------------------------------------------
490 subroutine bond_regular
494 include 'COMMON.LOCAL'
495 include 'COMMON.INTERACT'
496 include 'COMMON.CHAIN'
501 vbld(i+1+nres)=dsc(iabs(itype(i+1)))
502 vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
503 c print *,vbld(i+1),vbld(i+1+nres)
505 c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
511 vbld_inv(i1)=vbld_inv(i1)*2
514 vbld(i2+1)=vbld(i2+1)/2
515 vbld_inv(i2+1)=vbld_inv(i2+1)*2
520 c---------------------------------------------------------------------------
521 subroutine readpdb_template(k)
522 C Read the PDB file for read_constr_homology with read2sigma
523 C and convert the peptide geometry into virtual-chain geometry.
524 implicit real*8 (a-h,o-z)
526 include 'DIMENSIONS.ZSCOPT'
527 include 'COMMON.LOCAL'
529 include 'COMMON.CHAIN'
530 include 'COMMON.INTERACT'
531 include 'COMMON.IOUNITS'
533 include 'COMMON.NAMES'
534 include 'COMMON.CONTROL'
535 include 'COMMON.SETUP'
536 integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
537 logical lprn /.false./,fail
538 double precision e1(3),e2(3),e3(3)
539 double precision dcj,efree_temp
543 double precision sccor(3,20)
544 integer rescode,iterter(maxres)
551 c write (2,*) "UNRES_PDB",unres_pdb
559 read (ipdbin,'(a80)',end=10) card
560 if (card(:3).eq.'END') then
562 else if (card(:3).eq.'TER') then
565 itype(ires_old-1)=ntyp1
566 iterter(ires_old-1)=1
567 itype(ires_old)=ntyp1
570 c write (iout,*) "Chain ended",ires,ishift,ires_old
573 dc(j,ires)=sccor(j,iii)
576 call sccenter(ires,iii,sccor)
579 C Fish out the ATOM cards.
580 if (index(card(1:4),'ATOM').gt.0) then
581 read (card(12:16),*) atom
582 c write (iout,*) "! ",atom," !",ires
583 c if (atom.eq.'CA' .or. atom.eq.'CH3') then
584 read (card(23:26),*) ires
585 read (card(18:20),'(a3)') res
586 c write (iout,*) "ires",ires,ires-ishift+ishift1,
587 c & " ires_old",ires_old
588 c write (iout,*) "ishift",ishift," ishift1",ishift1
589 c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
590 if (ires-ishift+ishift1.ne.ires_old) then
591 C Calculate the CM of the preceding residue.
595 dc(j,ires)=sccor(j,iii)
598 call sccenter(ires_old,iii,sccor)
603 if (res.eq.'Cl-' .or. res.eq.'Na+') then
606 else if (ibeg.eq.1) then
607 c write (iout,*) "BEG ires",ires
609 if (res.ne.'GLY' .and. res.ne. 'ACE') then
613 ires=ires-ishift+ishift1
615 c write (iout,*) "ishift",ishift," ires",ires,
616 c & " ires_old",ires_old
617 c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
619 else if (ibeg.eq.2) then
621 ishift=-ires_old+ires-1
623 c write (iout,*) "New chain started",ires,ishift
626 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
627 ires=ires-ishift+ishift1
630 if (res.eq.'ACE' .or. res.eq.'NHE') then
633 itype(ires)=rescode(ires,res,0)
636 ires=ires-ishift+ishift1
638 c write (iout,*) "ires_old",ires_old," ires",ires
639 c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
642 c write (2,*) "ires",ires," res ",res," ity",ity
643 if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
644 & res.eq.'NHE'.and.atom(:2).eq.'HN') then
645 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
646 c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
648 write (iout,'(2i3,2x,a,3f8.3)')
649 & ires,itype(ires),res,(c(j,ires),j=1,3)
653 sccor(j,iii)=c(j,ires)
655 if (ishift.ne.0) then
656 ires_ca=ires+ishift-ishift1
660 c write (*,*) card(23:27),ires,itype(ires)
661 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
662 & atom.ne.'N' .and. atom.ne.'C' .and.
663 & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
664 & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
665 c write (iout,*) "sidechain ",atom
667 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
671 10 write (iout,'(a,i5)') ' Nres: ',ires
672 C Calculate dummy residue coordinates inside the "chain" of a multichain
676 c write (iout,*) i,itype(i),itype(i+1)
677 if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
678 if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
679 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
680 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
681 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
683 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
684 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
691 c(j,i)=c(j,i-1)-1.9d0*e2(j)
695 dcj=(c(j,i-2)-c(j,i-3))/2.0
696 if (dcj.eq.0) dcj=1.23591524223
701 else !itype(i+1).eq.ntyp1
703 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
704 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
711 c(j,i)=c(j,i+1)-1.9d0*e2(j)
715 dcj=(c(j,i+3)-c(j,i+2))/2.0
716 if (dcj.eq.0) dcj=1.23591524223
721 endif !itype(i+1).eq.ntyp1
722 endif !itype.eq.ntyp1
724 C Calculate the CM of the last side chain.
727 dc(j,ires)=sccor(j,iii)
730 call sccenter(ires,iii,sccor)
734 if (itype(nres).ne.10) then
738 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
739 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
746 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
750 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
751 if (dcj.eq.0) dcj=1.23591524223
752 c(j,nres)=c(j,nres-1)+dcj
753 c(j,2*nres)=c(j,nres)
764 c(j,2*nres)=c(j,nres)
766 if (itype(1).eq.ntyp1) then
770 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
771 call refsys(2,3,4,e1,e2,e3,fail)
778 c(j,1)=c(j,2)-1.9d0*e2(j)
782 dcj=(c(j,4)-c(j,3))/2.0
788 C Copy the coordinates to reference coordinates
794 C Calculate internal coordinates.
795 if (out_template_coord) then
797 & "Cartesian coordinates of the reference structure"
798 write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
799 & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
801 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
802 & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
803 & (c(j,ires+nres),j=1,3)
806 C Calculate internal coordinates.
807 c call int_from_cart1(.false.)
808 call int_from_cart(.true.,.true.)
809 call sc_loc_geom(.true.)
816 dc(j,i)=c(j,i+1)-c(j,i)
817 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
822 dc(j,i+nres)=c(j,i+nres)-c(j,i)
823 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
825 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
831 cref(j,i+nres)=c(j,i+nres)