Merge branch 'devel' of mmka.chem.univ.gda.pl:unres into devel
[unres.git] / source / wham / src-M / 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)=21
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)=21
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.21) 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)=21
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.21) 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.21) 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       do j=1,3
183       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
184       chain_rep(j,chain_length+nres,symetr)
185      &=chain_rep(j,chain_length+nres,1)
186       enddo
187 c diagnostic
188
189 c diagnostic
190 c       write (iout,*) "spraw lancuchy",chain_length,symetr
191 c       do i=1,24
192 c         do kkk=1,chain_length
193 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
194 c         enddo
195 c        enddo
196 c enddiagnostic       
197 C makes copy of chains
198 c      write (iout,*) "symetr", symetr
199
200       if (symetr.gt.1) then
201        call permut(symetr)
202        nperm=1
203        do i=1,symetr
204        nperm=nperm*i
205        enddo
206 c       do i=1,nperm
207 c       write(iout,*) "tabperm", (tabperm(i,kkk),kkk=1,4)
208 c       enddo
209        do i=1,nperm
210         cou=0
211         do kkk=1,symetr
212          icha=tabperm(i,kkk)
213 c         write (iout,*) i,icha
214          do lll=1,chain_length
215           cou=cou+1
216            if (cou.le.nres) then
217            do j=1,3
218             kupa=mod(lll,chain_length)
219             iprzes=(kkk-1)*chain_length+lll
220             if (kupa.eq.0) kupa=chain_length
221 c            write (iout,*) "kupa", kupa
222             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
223             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
224           enddo
225           endif
226          enddo
227         enddo
228        enddo
229        endif
230 C-koniec robienia kopidm
231       do kkk=1,nperm
232       write (iout,*) "nowa struktura", nperm
233       do i=1,nres
234         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
235      &cref(2,i,kkk),
236      &cref(3,i,kkk),cref(1,nres+i,kkk),
237      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
238       enddo
239   100 format (//'              alpha-carbon coordinates       ',
240      &          '     centroid coordinates'/
241      1          '       ', 6X,'X',11X,'Y',11X,'Z',
242      &                          10X,'X',11X,'Y',11X,'Z')
243   110 format (a,'(',i3,')',6f12.5)
244
245       enddo
246
247       ishift_pdb=ishift
248       return
249       end
250 c---------------------------------------------------------------------------
251       subroutine int_from_cart(lside,lprn)
252       implicit none
253       include 'DIMENSIONS'
254       include 'DIMENSIONS.ZSCOPT'
255       include 'COMMON.LOCAL'
256       include 'COMMON.VAR'
257       include 'COMMON.CHAIN'
258       include 'COMMON.INTERACT'
259       include 'COMMON.IOUNITS'
260       include 'COMMON.GEO'
261       include 'COMMON.NAMES'
262       character*3 seq,atom,res
263       character*80 card
264       double precision sccor(3,20)
265       integer rescode
266       double precision dist,alpha,beta,di
267       integer i,j,iti
268       logical lside,lprn
269       if (lprn) then 
270         write (iout,'(/a)') 
271      &  'Internal coordinates calculated from crystal structure.'
272         if (lside) then 
273           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
274      & '       Phi','    Dsc_id','       Dsc','     Alpha',
275      & '     Omega'
276         else 
277           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
278      & '       Phi'
279         endif
280       endif
281       do i=2,nres
282         iti=itype(i)
283         write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
284         if (itype(i-1).ne.21 .and. itype(i).ne.21 .and.
285      &    (dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0)) then
286           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
287           stop
288         endif
289         theta(i+1)=alpha(i-1,i,i+1)
290         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
291       enddo
292       if (itype(1).eq.21) then
293         do j=1,3
294           c(j,1)=c(j,2)+(c(j,3)-c(j,4))
295         enddo
296       endif
297       if (itype(nres).eq.21) then
298         do j=1,3
299           c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
300         enddo
301       endif
302       if (lside) then
303         do i=2,nres-1
304           do j=1,3
305             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
306           enddo
307           iti=itype(i)
308           di=dist(i,nres+i)
309           if (iti.ne.10) then
310             alph(i)=alpha(nres+i,i,maxres2)
311             omeg(i)=beta(nres+i,i,maxres2,i+1)
312           endif
313           if (lprn)
314      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
315      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
316      &    rad2deg*alph(i),rad2deg*omeg(i)
317         enddo
318       else if (lprn) then
319         do i=2,nres
320           iti=itype(i)
321           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
322      &    rad2deg*theta(i),rad2deg*phi(i)
323         enddo
324       endif
325       return
326       end
327 c---------------------------------------------------------------------------
328       subroutine sccenter(ires,nscat,sccor)
329       implicit none
330       include 'DIMENSIONS'
331       include 'COMMON.CHAIN'
332       integer ires,nscat,i,j
333       double precision sccor(3,20),sccmj
334       do j=1,3
335         sccmj=0.0D0
336         do i=1,nscat
337           sccmj=sccmj+sccor(j,i) 
338         enddo
339         dc(j,ires)=sccmj/nscat
340       enddo
341       return
342       end