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