Merge branch 'homology' of mmka.chem.univ.gda.pl:unres into homology
[unres.git] / source / unres / src_MD-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 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       integer rescode
21       logical fail
22       ibeg=1
23       lsecondary=.false.
24       nhfrag=0
25       nbfrag=0
26       do
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') then
46           goto 10
47         else if (card(:3).eq.'TER') then
48 C End current chain
49           ires_old=ires+1 
50           itype(ires_old)=21
51           ibeg=2
52 c          write (iout,*) "Chain ended",ires,ishift,ires_old
53           if (unres_pdb) then
54             do j=1,3
55               dc(j,ires)=sccor(j,iii)
56             enddo
57           else 
58             call sccenter(ires,iii,sccor)
59           endif
60         endif
61 C Fish out the ATOM cards.
62         if (index(card(1:4),'ATOM').gt.0) then  
63           read (card(14:16),'(a3)') atom
64           if (atom.eq.'CA' .or. atom.eq.'CH3') then
65 C Calculate the CM of the preceding residue.
66             if (ibeg.eq.0) then
67               if (unres_pdb) then
68                 do j=1,3
69                   dc(j,ires+nres)=sccor(j,iii)
70                 enddo
71               else
72                 call sccenter(ires,iii,sccor)
73               endif
74             endif
75 C Start new residue.
76 c            write (iout,'(a80)') card
77             read (card(24:26),*) ires
78             read (card(18:20),'(a3)') res
79             if (ibeg.eq.1) then
80               ishift=ires-1
81               if (res.ne.'GLY' .and. res.ne. 'ACE') then
82                 ishift=ishift-1
83                 itype(1)=21
84               endif
85 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
86               ibeg=0          
87             else if (ibeg.eq.2) then
88 c Start a new chain
89               ishift=-ires_old+ires-1
90 c              write (iout,*) "New chain started",ires,ishift
91               ibeg=0
92             endif
93             ires=ires-ishift
94 c            write (2,*) "ires",ires," ishift",ishift
95             if (res.eq.'ACE') then
96               ity=10
97             else
98               itype(ires)=rescode(ires,res,0)
99             endif
100             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
101             if(me.eq.king.or..not.out1file)
102      &       write (iout,'(2i3,2x,a,3f8.3)') 
103      &       ires,itype(ires),res,(c(j,ires),j=1,3)
104             iii=1
105             do j=1,3
106               sccor(j,iii)=c(j,ires)
107             enddo
108           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and. 
109      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
110             iii=iii+1
111             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
112           endif
113         endif
114       enddo
115    10 if(me.eq.king.or..not.out1file) 
116      & write (iout,'(a,i5)') ' Nres: ',ires
117 C Calculate dummy residue coordinates inside the "chain" of a multichain
118 C system
119       nres=ires
120       do i=2,nres-1
121 c        write (iout,*) i,itype(i)
122         if (itype(i).eq.21) then
123 c          write (iout,*) "dummy",i,itype(i)
124           do j=1,3
125             c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
126 c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
127             dc(j,i)=c(j,i)
128           enddo
129         endif
130       enddo
131 C Calculate the CM of the last side chain.
132       if (unres_pdb) then
133         do j=1,3
134           dc(j,ires)=sccor(j,iii)
135         enddo
136       else 
137         call sccenter(ires,iii,sccor)
138       endif
139       nsup=nres
140       nstart_sup=1
141       if (itype(nres).ne.10) then
142         nres=nres+1
143         itype(nres)=21
144         if (unres_pdb) then
145 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
146           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
147           if (fail) then
148             e2(1)=0.0d0
149             e2(2)=1.0d0
150             e2(3)=0.0d0
151           endif
152           do j=1,3
153             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
154           enddo
155         else
156         do j=1,3
157           dcj=c(j,nres-2)-c(j,nres-3)
158           c(j,nres)=c(j,nres-1)+dcj
159           c(j,2*nres)=c(j,nres)
160         enddo
161         endif
162       endif
163       do i=2,nres-1
164         do j=1,3
165           c(j,i+nres)=dc(j,i)
166         enddo
167       enddo
168       do j=1,3
169         c(j,nres+1)=c(j,1)
170         c(j,2*nres)=c(j,nres)
171       enddo
172       if (itype(1).eq.21) then 
173         nsup=nsup-1
174         nstart_sup=2
175         if (unres_pdb) then
176 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
177           call refsys(2,3,4,e1,e2,e3,fail)
178           if (fail) then
179             e2(1)=0.0d0
180             e2(2)=1.0d0
181             e2(3)=0.0d0
182           endif
183           do j=1,3
184             c(j,1)=c(j,2)-3.8d0*e2(j)
185           enddo
186         else
187         do j=1,3
188           dcj=c(j,4)-c(j,3)
189           c(j,1)=c(j,2)-dcj
190           c(j,nres+1)=c(j,1)
191         enddo
192         endif
193       endif
194 C Calculate internal coordinates.
195       if(me.eq.king.or..not.out1file)then
196        do ires=1,nres
197         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') 
198      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
199      &    (c(j,nres+ires),j=1,3)
200        enddo
201       endif
202       call int_from_cart(.true.,.false.)
203       call sc_loc_geom(.true.)
204       do i=1,nres
205         thetaref(i)=theta(i)
206         phiref(i)=phi(i)
207       enddo
208       do i=1,nres-1
209         do j=1,3
210           dc(j,i)=c(j,i+1)-c(j,i)
211           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
212         enddo
213       enddo
214       do i=2,nres-1
215         do j=1,3
216           dc(j,i+nres)=c(j,i+nres)-c(j,i)
217           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
218         enddo
219 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
220 c     &   vbld_inv(i+nres)
221       enddo
222 c      call chainbuild
223 C Copy the coordinates to reference coordinates
224 C Splits to single chain if occurs
225       kkk=1
226       lll=0
227       cou=1
228       do i=1,nres
229       lll=lll+1
230 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
231       if (i.gt.1) then
232       if ((itype(i-1).eq.21)) then
233       chain_length=lll-1
234       kkk=kkk+1
235 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
236       lll=1
237       endif
238       endif
239         do j=1,3
240           cref(j,i,cou)=c(j,i)
241           cref(j,i+nres,cou)=c(j,i+nres)
242           if (i.le.nres) then
243           chain_rep(j,lll,kkk)=c(j,i)
244           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
245           endif
246          enddo
247       enddo
248       do j=1,3
249       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
250       chain_rep(j,chain_length+nres,symetr)
251      &=chain_rep(j,chain_length+nres,1)
252       enddo
253 c diagnostic
254 c       write (iout,*) "spraw lancuchy",chain_length,symetr
255 c       do i=1,4
256 c         do kkk=1,chain_length
257 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
258 c         enddo
259 c        enddo
260 c enddiagnostic       
261 C makes copy of chains
262         write (iout,*) "symetr", symetr
263        
264       if (symetr.gt.1) then
265        call permut(symetr)
266        nperm=1
267        do i=1,symetr
268        nperm=nperm*i
269        enddo
270        do i=1,nperm
271        write(iout,*) (tabperm(i,kkk),kkk=1,4)
272        enddo
273        do i=1,nperm
274         cou=0
275         do kkk=1,symetr
276          icha=tabperm(i,kkk)
277 c         write (iout,*) i,icha
278          do lll=1,chain_length
279           cou=cou+1
280            if (cou.le.nres) then
281            do j=1,3
282             kupa=mod(lll,chain_length)
283             iprzes=(kkk-1)*chain_length+lll
284             if (kupa.eq.0) kupa=chain_length
285 c            write (iout,*) "kupa", kupa
286             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
287             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
288           enddo
289           endif
290          enddo
291         enddo
292        enddo
293        endif
294 C-koniec robienia kopii
295 c diag
296       do kkk=1,nperm
297       write (iout,*) "nowa struktura", nperm
298       do i=1,nres
299         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
300      &cref(2,i,kkk),
301      &cref(3,i,kkk),cref(1,nres+i,kkk),
302      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
303       enddo
304   100 format (//'              alpha-carbon coordinates       ',
305      &          '     centroid coordinates'/
306      1          '       ', 6X,'X',11X,'Y',11X,'Z',
307      &                          10X,'X',11X,'Y',11X,'Z')
308   110 format (a,'(',i3,')',6f12.5)
309       
310       enddo
311 cc enddiag
312       do j=1,nbfrag     
313         do i=1,4                                                       
314          bfrag(i,j)=bfrag(i,j)-ishift
315         enddo
316       enddo
317
318       do j=1,nhfrag
319         do i=1,2
320          hfrag(i,j)=hfrag(i,j)-ishift
321         enddo
322       enddo
323
324       return
325       end
326 c---------------------------------------------------------------------------
327       subroutine int_from_cart(lside,lprn)
328       implicit real*8 (a-h,o-z)
329       include 'DIMENSIONS'
330 #ifdef MPI
331       include "mpif.h"
332 #endif 
333       include 'COMMON.LOCAL'
334       include 'COMMON.VAR'
335       include 'COMMON.CHAIN'
336       include 'COMMON.INTERACT'
337       include 'COMMON.IOUNITS'
338       include 'COMMON.GEO'
339       include 'COMMON.NAMES'
340       include 'COMMON.CONTROL'
341       include 'COMMON.SETUP'
342       character*3 seq,atom,res
343       character*80 card
344       dimension sccor(3,20)
345       integer rescode
346       logical lside,lprn
347 #ifdef MPI
348       if(me.eq.king.or..not.out1file)then
349 #endif
350        if (lprn) then 
351         write (iout,'(/a)') 
352      &  'Internal coordinates calculated from crystal structure.'
353         if (lside) then 
354           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
355      & '       Phi','    Dsc_id','       Dsc','     Alpha',
356      & '     Omega'
357         else 
358           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
359      & '       Phi'
360         endif
361        endif
362 #ifdef MPI
363       endif
364 #endif
365       do i=1,nres-1
366         iti=itype(i)
367         if (iti.ne.21 .and. itype(i+1).ne.21 .and. 
368      &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
369           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
370 ctest          stop
371         endif
372         vbld(i+1)=dist(i,i+1)
373         vbld_inv(i+1)=1.0d0/vbld(i+1)
374         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
375         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
376       enddo
377 c      if (unres_pdb) then
378 c        if (itype(1).eq.21) then
379 c          theta(3)=90.0d0*deg2rad
380 c          phi(4)=180.0d0*deg2rad
381 c          vbld(2)=3.8d0
382 c          vbld_inv(2)=1.0d0/vbld(2)
383 c        endif
384 c        if (itype(nres).eq.21) then
385 c          theta(nres)=90.0d0*deg2rad
386 c          phi(nres)=180.0d0*deg2rad
387 c          vbld(nres)=3.8d0
388 c          vbld_inv(nres)=1.0d0/vbld(2)
389 c        endif
390 c      endif
391       if (lside) then
392         do i=2,nres-1
393           do j=1,3
394             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
395      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
396           enddo
397           iti=itype(i)
398           di=dist(i,nres+i)
399           vbld(i+nres)=di
400           if (itype(i).ne.10) then
401             vbld_inv(i+nres)=1.0d0/di
402           else
403             vbld_inv(i+nres)=0.0d0
404           endif
405           if (iti.ne.10) then
406             alph(i)=alpha(nres+i,i,maxres2)
407             omeg(i)=beta(nres+i,i,maxres2,i+1)
408           endif
409           if(me.eq.king.or..not.out1file)then
410            if (lprn)
411      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
412      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
413      &     rad2deg*alph(i),rad2deg*omeg(i)
414           endif
415         enddo
416       else if (lprn) then
417         do i=2,nres
418           iti=itype(i)
419           if(me.eq.king.or..not.out1file)
420      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
421      &     rad2deg*theta(i),rad2deg*phi(i)
422         enddo
423       endif
424       return
425       end
426 c-------------------------------------------------------------------------------
427       subroutine sc_loc_geom(lprn)
428       implicit real*8 (a-h,o-z)
429       include 'DIMENSIONS'
430 #ifdef MPI
431       include "mpif.h"
432 #endif 
433       include 'COMMON.LOCAL'
434       include 'COMMON.VAR'
435       include 'COMMON.CHAIN'
436       include 'COMMON.INTERACT'
437       include 'COMMON.IOUNITS'
438       include 'COMMON.GEO'
439       include 'COMMON.NAMES'
440       include 'COMMON.CONTROL'
441       include 'COMMON.SETUP'
442       double precision x_prime(3),y_prime(3),z_prime(3)
443       logical lprn
444       do i=1,nres-1
445         do j=1,3
446           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
447         enddo
448       enddo
449       do i=2,nres-1
450         if (itype(i).ne.10 .and. itype(i).ne.21) then
451           do j=1,3
452             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
453           enddo
454         else
455           do j=1,3
456             dc_norm(j,i+nres)=0.0d0
457           enddo
458         endif
459       enddo
460       do i=2,nres-1
461         costtab(i+1) =dcos(theta(i+1))
462         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
463         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
464         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
465         cosfac2=0.5d0/(1.0d0+costtab(i+1))
466         cosfac=dsqrt(cosfac2)
467         sinfac2=0.5d0/(1.0d0-costtab(i+1))
468         sinfac=dsqrt(sinfac2)
469         it=itype(i)
470         if (it.ne.10 .and. itype(i).ne.21) then
471 c
472 C  Compute the axes of tghe local cartesian coordinates system; store in
473 c   x_prime, y_prime and z_prime 
474 c
475         do j=1,3
476           x_prime(j) = 0.00
477           y_prime(j) = 0.00
478           z_prime(j) = 0.00
479         enddo
480         do j = 1,3
481           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
482           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
483         enddo
484         call vecpr(x_prime,y_prime,z_prime)
485 c
486 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
487 C to local coordinate system. Store in xx, yy, zz.
488 c
489         xx=0.0d0
490         yy=0.0d0
491         zz=0.0d0
492         do j = 1,3
493           xx = xx + x_prime(j)*dc_norm(j,i+nres)
494           yy = yy + y_prime(j)*dc_norm(j,i+nres)
495           zz = zz + z_prime(j)*dc_norm(j,i+nres)
496         enddo
497
498         xxref(i)=xx
499         yyref(i)=yy
500         zzref(i)=zz
501         else
502         xxref(i)=0.0d0
503         yyref(i)=0.0d0
504         zzref(i)=0.0d0
505         endif
506       enddo
507       if (lprn) then
508         do i=2,nres
509           iti=itype(i)
510 #ifdef MPI
511           if(me.eq.king.or..not.out1file)
512      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
513      &      yyref(i),zzref(i)
514 #else
515           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
516      &     zzref(i)
517 #endif
518         enddo
519       endif
520       return
521       end
522 c---------------------------------------------------------------------------
523       subroutine sccenter(ires,nscat,sccor)
524       implicit real*8 (a-h,o-z)
525       include 'DIMENSIONS'
526       include 'COMMON.CHAIN'
527       dimension sccor(3,20)
528       do j=1,3
529         sccmj=0.0D0
530         do i=1,nscat
531           sccmj=sccmj+sccor(j,i) 
532         enddo
533         dc(j,ires)=sccmj/nscat
534       enddo
535       return
536       end
537 c---------------------------------------------------------------------------
538       subroutine bond_regular
539       implicit real*8 (a-h,o-z)
540       include 'DIMENSIONS'   
541       include 'COMMON.VAR'
542       include 'COMMON.LOCAL'      
543       include 'COMMON.CALC'
544       include 'COMMON.INTERACT'
545       include 'COMMON.CHAIN'
546       do i=1,nres-1
547        vbld(i+1)=vbl
548        vbld_inv(i+1)=1.0d0/vbld(i+1)
549        vbld(i+1+nres)=dsc(itype(i+1))
550        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
551 c       print *,vbld(i+1),vbld(i+1+nres)
552       enddo
553       return
554       end
555