725065a67435df41a83f5ed93be88a61da0b329f
[unres.git] / readpdb.f
1       subroutine readpdb
2 C Read the PDB file and convert the peptide geometry into virtual-chain 
3 C geometry.
4       implicit real*8 (a-h,o-z)
5       include 'DIMENSIONS'
6       include 'COMMON.LOCAL'
7       include 'COMMON.VAR'
8       include 'COMMON.CHAIN'
9       include 'COMMON.INTERACT'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.GEO'
12       include 'COMMON.NAMES'
13       character*3 seq,atom,res
14       character*80 card
15       dimension sccor(3,20)
16       integer rescode
17       call permut(symetr)
18       ibeg=1
19       do
20         read (ipdbin,'(a80)',end=10) card
21         if (card(:3).eq.'END') then
22           goto 10
23         else if (card(:3).eq.'TER') then
24 C End current chain
25           ires_old=ires+1 
26           itype(ires_old)=ntyp1
27           ibeg=2
28 c          write (iout,*) "Chain ended",ires,ishift,ires_old
29           call sccenter(ires,iii,sccor)
30         endif
31 C Fish out the ATOM cards.
32         if (index(card(1:4),'ATOM').gt.0) then  
33           read (card(14:16),'(a3)') atom
34           if (atom.eq.'CA' .or. atom.eq.'CH3') then
35 C Calculate the CM of the preceding residue.
36             if (ibeg.eq.0) then
37               call sccenter(ires,iii,sccor)
38             endif
39 C Start new residue.
40 c            write (iout,'(a80)') card
41             read (card(24:26),*) ires
42             read (card(18:20),'(a3)') res
43             if (ibeg.eq.1) then
44               ishift=ires-1
45               if (res.ne.'GLY' .and. res.ne. 'ACE') then
46                 ishift=ishift-1
47                 itype(1)=ntyp1
48               endif
49 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
50               ibeg=0          
51             else if (ibeg.eq.2) then
52 c Start a new chain
53               ishift=-ires_old+ires-1
54 c              write (iout,*) "New chain started",ires,ishift
55               ibeg=0
56             endif
57             ires=ires-ishift
58 c            write (2,*) "ires",ires," ishift",ishift
59             if (res.eq.'ACE') then
60               ity=10
61             else
62               itype(ires)=rescode(ires,res,0)
63             endif
64             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
65             write (iout,'(2i3,2x,a,3f8.3)') 
66      &       ires,itype(ires),res,(c(j,ires),j=1,3)
67             iii=1
68             do j=1,3
69               sccor(j,iii)=c(j,ires)
70             enddo
71           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and. 
72      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
73             iii=iii+1
74             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
75           endif
76         endif
77       enddo
78    10 write (iout,'(a,i5)') ' Nres: ',ires
79 C Calculate dummy residue coordinates inside the "chain" of a multichain
80 C system
81       nres=ires
82       do i=2,nres-1
83 c        write (iout,*) i,itype(i)
84         if (itype(i).eq.ntyp1) then
85 c          write (iout,*) "dummy",i,itype(i)
86           do j=1,3
87             c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
88 c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
89             dc(j,i)=c(j,i)
90           enddo
91         endif
92       enddo
93 C Calculate the CM of the last side chain.
94       call sccenter(ires,iii,sccor)
95       nsup=nres
96       nstart_sup=1
97       if (itype(nres).ne.10) then
98         nres=nres+1
99         itype(nres)=ntyp1
100         do j=1,3
101           dcj=c(j,nres-2)-c(j,nres-3)
102           c(j,nres)=c(j,nres-1)+dcj
103           c(j,2*nres)=c(j,nres)
104         enddo
105       endif
106       do i=2,nres-1
107         do j=1,3
108           c(j,i+nres)=dc(j,i)
109         enddo
110       enddo
111       do j=1,3
112         c(j,nres+1)=c(j,1)
113         c(j,2*nres)=c(j,nres)
114       enddo
115       if (itype(1).eq.ntyp1) then
116         nsup=nsup-1
117         nstart_sup=2
118         do j=1,3
119           dcj=c(j,4)-c(j,3)
120           c(j,1)=c(j,2)-dcj
121           c(j,nres+1)=c(j,1)
122         enddo
123       endif
124 C Calculate internal coordinates.
125       do ires=1,nres
126         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') 
127      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
128      &    (c(j,nres+ires),j=1,3)
129       enddo
130       call int_from_cart(.true.,.false.)
131       do i=1,nres-1
132         do j=1,3
133           dc(j,i)=c(j,i+1)-c(j,i)
134           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
135         enddo
136       enddo
137       do i=2,nres-1
138         do j=1,3
139           dc(j,i+nres)=c(j,i+nres)-c(j,i)
140           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
141         enddo
142 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
143 c     &   vbld_inv(i+nres)
144       enddo
145 c      call chainbuild
146 C Copy the coordinates to reference coordinates
147 c      do i=1,2*nres
148 c        do j=1,3
149 c          cref(j,i)=c(j,i)
150 c        enddo
151 c      enddo
152
153       kkk=1
154       lll=0
155       cou=1
156       do i=1,nres
157       lll=lll+1
158 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
159       if (i.gt.1) then
160       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
161       chain_length=lll-1
162       kkk=kkk+1
163 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
164       lll=1
165       endif
166       endif
167         do j=1,3
168           cref(j,i,cou)=c(j,i)
169           cref(j,i+nres,cou)=c(j,i+nres)
170           if (i.le.nres) then
171           chain_rep(j,lll,kkk)=c(j,i)
172           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
173           endif
174          enddo
175       enddo
176       do j=1,3
177       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
178       chain_rep(j,chain_length+nres,symetr)
179      &=chain_rep(j,chain_length+nres,1)
180       enddo
181
182       if (symetr.gt.1) then
183        call permut(symetr)
184        nperm=1
185        do i=1,symetr
186        nperm=nperm*i
187        enddo
188 c       do i=1,nperm
189 c       write(iout,*) "tabperm", (tabperm(i,kkk),kkk=1,4)
190 c       enddo
191        do i=1,nperm
192         cou=0
193         do kkk=1,symetr
194          icha=tabperm(i,kkk)
195 c         write (iout,*) i,icha
196          do lll=1,chain_length
197           cou=cou+1
198            if (cou.le.nres) then
199            do j=1,3
200             kupa=mod(lll,chain_length)
201             iprzes=(kkk-1)*chain_length+lll
202             if (kupa.eq.0) kupa=chain_length
203 c            write (iout,*) "kupa", kupa
204             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
205             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
206           enddo
207           endif
208          enddo
209         enddo
210        enddo
211        endif
212
213 C-koniec robienia kopidm
214       do kkk=1,nperm
215       write (iout,*) "nowa struktura", nperm
216       do i=1,nres
217         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
218      &cref(2,i,kkk),
219      &cref(3,i,kkk),cref(1,nres+i,kkk),
220      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
221       enddo
222   100 format (//'              alpha-carbon coordinates       ',
223      &          '     centroid coordinates'/
224      1          '       ', 6X,'X',11X,'Y',11X,'Z',
225      &                          10X,'X',11X,'Y',11X,'Z')
226   110 format (a,'(',i3,')',6f12.5)
227        enddo
228
229
230       ishift_pdb=ishift
231       return
232       end
233 c---------------------------------------------------------------------------
234       subroutine int_from_cart(lside,lprn)
235       implicit real*8 (a-h,o-z)
236       include 'DIMENSIONS'
237       include 'COMMON.LOCAL'
238       include 'COMMON.VAR'
239       include 'COMMON.CHAIN'
240       include 'COMMON.INTERACT'
241       include 'COMMON.IOUNITS'
242       include 'COMMON.GEO'
243       include 'COMMON.NAMES'
244       character*3 seq,atom,res
245       character*80 card
246       dimension sccor(3,20)
247       integer rescode
248       logical lside,lprn
249       if (lprn) then 
250         write (iout,'(/a)') 
251      &  'Internal coordinates calculated from crystal structure.'
252         if (lside) then 
253           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
254      & '       Phi','    Dsc_id','       Dsc','     Alpha',
255      & '     Omega'
256         else 
257           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
258      & '       Phi'
259         endif
260       endif
261       call flush(iout)
262       do i=nnt+1,nct
263         iti=itype(i)
264 c        write (iout,*) i,dist(i,i-1)
265         if (dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0) then
266           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
267           stop
268         endif
269         theta(i+1)=alpha(i-1,i,i+1)
270         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
271       enddo
272       if (lside) then
273         do i=2,nres-1
274           do j=1,3
275             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
276           enddo
277           iti=itype(i)
278           di=dist(i,nres+i)
279           if (iti.ne.10) then
280             alph(i)=alpha(nres+i,i,maxres2)
281             omeg(i)=beta(nres+i,i,maxres2,i+1)
282           endif
283           if (lprn)
284      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
285      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,rad2deg*alph(i),
286      &    rad2deg*omeg(i)
287         enddo
288       else if (lprn) then
289         do i=2,nres
290           iti=itype(i)
291           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
292      &    rad2deg*theta(i),rad2deg*phi(i)
293         enddo
294       endif
295       return
296       end
297 c---------------------------------------------------------------------------
298       subroutine sccenter(ires,nscat,sccor)
299       implicit real*8 (a-h,o-z)
300       include 'DIMENSIONS'
301       include 'COMMON.CHAIN'
302       dimension sccor(3,20)
303       do j=1,3
304         sccmj=0.0D0
305         do i=1,nscat
306           sccmj=sccmj+sccor(j,i) 
307         enddo
308         dc(j,ires)=sccmj/nscat
309       enddo
310       return
311       end