Fixed the addition of dummy chain ends to structures read from UNRES PDB files
[unres.git] / source / unres / src_CSA / 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       include 'COMMON.CONTROL'
14       include 'COMMON.DISTFIT'
15       include 'COMMON.SETUP'
16       character*3 seq,atom,res
17       character*80 card
18       dimension sccor(3,20)
19       double precision e1(3),e2(3),e3(3)
20       logical fail
21       integer rescode
22       ibeg=1
23       lsecondary=.false.
24       nhfrag=0
25       nbfrag=0
26       do i=1,10000
27         read (ipdbin,'(a80)',end=10) card
28         if (card(:5).eq.'HELIX') then
29          nhfrag=nhfrag+1
30          lsecondary=.true.
31          read(card(22:25),*) hfrag(1,nhfrag)
32          read(card(34:37),*) hfrag(2,nhfrag)
33         endif
34         if (card(:5).eq.'SHEET') then
35          nbfrag=nbfrag+1
36          lsecondary=.true.
37          read(card(24:26),*) bfrag(1,nbfrag)
38          read(card(35:37),*) bfrag(2,nbfrag)
39 crc----------------------------------------
40 crc  to be corrected !!!
41          bfrag(3,nbfrag)=bfrag(1,nbfrag)
42          bfrag(4,nbfrag)=bfrag(2,nbfrag)
43 crc----------------------------------------
44         endif
45         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
46 C Fish out the ATOM cards.
47         if (index(card(1:4),'ATOM').gt.0) then  
48           read (card(14:16),'(a3)') atom
49           if (atom.eq.'CA' .or. atom.eq.'CH3') then
50 C Calculate the CM of the preceding residue.
51             if (ibeg.eq.0) then
52               if (unres_pdb) then
53                 do j=1,3
54                   dc(j,ires+nres)=sccor(j,iii)
55                 enddo
56               else
57                 call sccenter(ires,iii,sccor)
58               endif
59             endif
60 C Start new residue.
61             read (card(24:26),*) ires
62             read (card(18:20),'(a3)') res
63             if (ibeg.eq.1) then
64               ishift=ires-1
65               if (res.ne.'GLY' .and. res.ne. 'ACE') then
66                 ishift=ishift-1
67                 itype(1)=21
68               endif
69               ibeg=0          
70             endif
71             ires=ires-ishift
72             if (res.eq.'ACE') then
73               ity=10
74             else
75               itype(ires)=rescode(ires,res,0)
76             endif
77             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
78 c            if(me.eq.king.or..not.out1file)
79 c     &       write (iout,'(2i3,2x,a,3f8.3)') 
80 c     &       ires,itype(ires),res,(c(j,ires),j=1,3)
81             iii=1
82             do j=1,3
83               sccor(j,iii)=c(j,ires)
84             enddo
85           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and. 
86      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
87             iii=iii+1
88             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
89           endif
90         endif
91       enddo
92    10 if(me.eq.king.or..not.out1file) 
93      & write (iout,'(a,i5)') ' Nres: ',ires
94 C Calculate the CM of the last side chain.
95       if (unres_pdb) then
96         do j=1,3
97           dc(j,ires+nres)=sccor(j,iii)
98         enddo
99       else 
100         call sccenter(ires,iii,sccor)
101       endif
102       nres=ires
103       nsup=nres
104       nstart_sup=1
105       if (itype(nres).ne.10) then
106         nres=nres+1
107         itype(nres)=21
108         if (unres_pdb) then
109 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
110           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
111           if (fail) then
112             e2(1)=0.0d0
113             e2(2)=1.0d0
114             e2(3)=0.0d0
115           endif
116           do j=1,3
117             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
118           enddo
119         else
120         do j=1,3
121           dcj=c(j,nres-2)-c(j,nres-3)
122           c(j,nres)=c(j,nres-1)+dcj
123           c(j,2*nres)=c(j,nres)
124         enddo
125         endif
126       endif
127       do i=2,nres-1
128         do j=1,3
129           c(j,i+nres)=dc(j,i)
130         enddo
131       enddo
132       do j=1,3
133         c(j,nres+1)=c(j,1)
134         c(j,2*nres)=c(j,nres)
135       enddo
136       if (itype(1).eq.21) then
137         nsup=nsup-1
138         nstart_sup=2
139         if (unres_pdb) then
140 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
141           call refsys(2,3,4,e1,e2,e3,fail)
142           if (fail) then
143             e2(1)=0.0d0
144             e2(2)=1.0d0
145             e2(3)=0.0d0
146           endif
147           do j=1,3
148             c(j,1)=c(j,2)-3.8d0*e2(j)
149           enddo
150         else
151         do j=1,3
152           dcj=c(j,4)-c(j,3)
153           c(j,1)=c(j,2)-dcj
154           c(j,nres+1)=c(j,1)
155         enddo
156         endif
157       endif
158 C Calculate internal coordinates.
159       if(me.eq.king.or..not.out1file)then
160        write (iout,'(a)') 
161      &   "Backbone and SC coordinates as read from the PDB"
162        do ires=1,nres
163         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') 
164      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
165      &    (c(j,nres+ires),j=1,3)
166        enddo
167       endif
168       call int_from_cart(.true.,.false.)
169       call sc_loc_geom(.false.)
170       do i=1,nres
171         thetaref(i)=theta(i)
172         phiref(i)=phi(i)
173       enddo
174       do i=1,nres-1
175         do j=1,3
176           dc(j,i)=c(j,i+1)-c(j,i)
177           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
178         enddo
179       enddo
180       do i=2,nres-1
181         do j=1,3
182           dc(j,i+nres)=c(j,i+nres)-c(j,i)
183           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
184         enddo
185 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
186 c     &   vbld_inv(i+nres)
187       enddo
188 c      call chainbuild
189 C Copy the coordinates to reference coordinates
190       do i=1,2*nres
191         do j=1,3
192           cref(j,i)=c(j,i)
193         enddo
194       enddo
195
196
197       do j=1,nbfrag     
198         do i=1,4                                                       
199          bfrag(i,j)=bfrag(i,j)-ishift
200         enddo
201       enddo
202
203       do j=1,nhfrag
204         do i=1,2
205          hfrag(i,j)=hfrag(i,j)-ishift
206         enddo
207       enddo
208
209       return
210       end
211 c---------------------------------------------------------------------------
212       subroutine int_from_cart(lside,lprn)
213       implicit real*8 (a-h,o-z)
214       include 'DIMENSIONS'
215 c      include "mpif.h"
216       include 'COMMON.LOCAL'
217       include 'COMMON.VAR'
218       include 'COMMON.CHAIN'
219       include 'COMMON.INTERACT'
220       include 'COMMON.IOUNITS'
221       include 'COMMON.GEO'
222       include 'COMMON.NAMES'
223       include 'COMMON.CONTROL'
224       include 'COMMON.SETUP'
225       character*3 seq,atom,res
226       character*80 card
227       dimension sccor(3,20)
228       integer rescode
229       logical lside,lprn
230       if(me.eq.king.or..not.out1file)then
231        if (lprn) then 
232         write (iout,'(/a)') 
233      &  'Internal coordinates calculated from crystal structure.'
234         if (lside) then 
235           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
236      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
237      & '     Beta '
238         else 
239           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
240      & '     Gamma'
241         endif
242        endif
243       endif
244       do i=1,nres-1
245         iti=itype(i)
246         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
247           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
248 ctest          stop
249         endif
250         vbld(i+1)=dist(i,i+1)
251         vbld_inv(i+1)=1.0d0/vbld(i+1)
252         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
253         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
254       enddo
255 c      if (unres_pdb) then
256 c        if (itype(1).eq.21) then
257 c          theta(3)=90.0d0*deg2rad
258 c          phi(4)=180.0d0*deg2rad
259 c          vbld(2)=3.8d0
260 c          vbld_inv(2)=1.0d0/vbld(2)
261 c        endif
262 c        if (itype(nres).eq.21) then
263 c          theta(nres)=90.0d0*deg2rad
264 c          phi(nres)=180.0d0*deg2rad
265 c          vbld(nres)=3.8d0
266 c          vbld_inv(nres)=1.0d0/vbld(2)
267 c        endif
268 c      endif
269       if (lside) then
270         do i=2,nres-1
271           do j=1,3
272             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
273      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
274           enddo
275           iti=itype(i)
276           di=dist(i,nres+i)
277 C 9/29/12 Adam: Correction for zero SC-SC bond length 
278           if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0) 
279      &     di=dsc(itype(i))
280           vbld(i+nres)=di
281           if (itype(i).ne.10) then
282             vbld_inv(i+nres)=1.0d0/di
283           else
284             vbld_inv(i+nres)=0.0d0
285           endif
286           if (iti.ne.10) then
287             alph(i)=alpha(nres+i,i,maxres2)
288             omeg(i)=beta(nres+i,i,maxres2,i+1)
289           endif
290           if(me.eq.king.or..not.out1file)then
291            if (lprn)
292      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
293      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
294      &     rad2deg*alph(i),rad2deg*omeg(i)
295           endif
296         enddo
297       else if (lprn) then
298         do i=2,nres
299           iti=itype(i)
300           if(me.eq.king.or..not.out1file)
301      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
302      &     rad2deg*theta(i),rad2deg*phi(i)
303         enddo
304       endif
305       return
306       end
307 c-------------------------------------------------------------------------------
308       subroutine sc_loc_geom(lprn)
309       implicit real*8 (a-h,o-z)
310       include 'DIMENSIONS'
311 c      include "mpif.h"
312       include 'COMMON.LOCAL'
313       include 'COMMON.VAR'
314       include 'COMMON.CHAIN'
315       include 'COMMON.INTERACT'
316       include 'COMMON.IOUNITS'
317       include 'COMMON.GEO'
318       include 'COMMON.NAMES'
319       include 'COMMON.CONTROL'
320       include 'COMMON.SETUP'
321       double precision x_prime(3),y_prime(3),z_prime(3)
322       logical lprn
323       do i=1,nres-1
324         do j=1,3
325           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
326         enddo
327       enddo
328       do i=2,nres-1
329         if (itype(i).ne.10) then
330           do j=1,3
331             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
332           enddo
333         else
334           do j=1,3
335             dc_norm(j,i+nres)=0.0d0
336           enddo
337         endif
338       enddo
339       do i=2,nres-1
340         costtab(i+1) =dcos(theta(i+1))
341         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
342         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
343         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
344         cosfac2=0.5d0/(1.0d0+costtab(i+1))
345         cosfac=dsqrt(cosfac2)
346         sinfac2=0.5d0/(1.0d0-costtab(i+1))
347         sinfac=dsqrt(sinfac2)
348         it=itype(i)
349         if (it.ne.10) then
350 c
351 C  Compute the axes of tghe local cartesian coordinates system; store in
352 c   x_prime, y_prime and z_prime 
353 c
354         do j=1,3
355           x_prime(j) = 0.00
356           y_prime(j) = 0.00
357           z_prime(j) = 0.00
358         enddo
359         do j = 1,3
360           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
361           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
362         enddo
363         call vecpr(x_prime,y_prime,z_prime)
364 c
365 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
366 C to local coordinate system. Store in xx, yy, zz.
367 c
368         xx=0.0d0
369         yy=0.0d0
370         zz=0.0d0
371         do j = 1,3
372           xx = xx + x_prime(j)*dc_norm(j,i+nres)
373           yy = yy + y_prime(j)*dc_norm(j,i+nres)
374           zz = zz + z_prime(j)*dc_norm(j,i+nres)
375         enddo
376
377         xxref(i)=xx
378         yyref(i)=yy
379         zzref(i)=zz
380         else
381         xxref(i)=0.0d0
382         yyref(i)=0.0d0
383         zzref(i)=0.0d0
384         endif
385       enddo
386       if (lprn) then
387         do i=2,nres
388           iti=itype(i)
389           if(me.eq.king.or..not.out1file)
390      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
391      &      yyref(i),zzref(i)
392         enddo
393       endif
394       return
395       end
396 c---------------------------------------------------------------------------
397       subroutine sccenter(ires,nscat,sccor)
398       implicit real*8 (a-h,o-z)
399       include 'DIMENSIONS'
400       include 'COMMON.CHAIN'
401       dimension sccor(3,20)
402       do j=1,3
403         sccmj=0.0D0
404         do i=1,nscat
405           sccmj=sccmj+sccor(j,i) 
406         enddo
407         dc(j,ires)=sccmj/nscat
408       enddo
409       return
410       end
411 c---------------------------------------------------------------------------
412       subroutine bond_regular
413       implicit real*8 (a-h,o-z)
414       include 'DIMENSIONS'   
415       include 'COMMON.VAR'
416       include 'COMMON.LOCAL'      
417       include 'COMMON.CALC'
418       include 'COMMON.INTERACT'
419       include 'COMMON.CHAIN'
420       do i=1,nres-1
421        vbld(i+1)=vbl
422        vbld_inv(i+1)=1.0d0/vbld(i+1)
423        vbld(i+1+nres)=dsc(itype(i+1))
424        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
425 c       print *,vbld(i+1),vbld(i+1+nres)
426       enddo
427       return
428       end