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