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