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