CMake makro cinfo.f fix
[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       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,nres
213       lll=lll+1
214 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
215       if (i.gt.1) then
216       if ((itype(i-1).eq.21)) then
217       chain_length=lll-1
218       kkk=kkk+1
219 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
220       lll=1
221       endif
222       endif
223         do j=1,3
224           cref(j,i,cou)=c(j,i)
225           cref(j,i+nres,cou)=c(j,i+nres)
226           if (i.le.nres) then
227           chain_rep(j,lll,kkk)=c(j,i)
228           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
229           endif
230          enddo
231       enddo
232       do j=1,3
233       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
234       chain_rep(j,chain_length+nres,symetr)
235      &=chain_rep(j,chain_length+nres,1)
236       enddo
237 c diagnostic
238 c       write (iout,*) "spraw lancuchy",chain_length,symetr
239 c       do i=1,4
240 c         do kkk=1,chain_length
241 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
242 c         enddo
243 c        enddo
244 c enddiagnostic       
245 C makes copy of chains
246         write (iout,*) "symetr", symetr
247        
248       if (symetr.gt.1) then
249        call permut(symetr)
250        nperm=1
251        do i=1,symetr
252        nperm=nperm*i
253        enddo
254        do i=1,nperm
255        write(iout,*) (tabperm(i,kkk),kkk=1,4)
256        enddo
257        do i=1,nperm
258         cou=0
259         do kkk=1,symetr
260          icha=tabperm(i,kkk)
261 c         write (iout,*) i,icha
262          do lll=1,chain_length
263           cou=cou+1
264            if (cou.le.nres) then
265            do j=1,3
266             kupa=mod(lll,chain_length)
267             iprzes=(kkk-1)*chain_length+lll
268             if (kupa.eq.0) kupa=chain_length
269 c            write (iout,*) "kupa", kupa
270             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
271             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
272           enddo
273           endif
274          enddo
275         enddo
276        enddo
277        endif
278 C-koniec robienia kopii
279 c diag
280       do kkk=1,nperm
281       write (iout,*) "nowa struktura", nperm
282       do i=1,nres
283         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
284      &cref(2,i,kkk),
285      &cref(3,i,kkk),cref(1,nres+i,kkk),
286      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
287       enddo
288   100 format (//'              alpha-carbon coordinates       ',
289      &          '     centroid coordinates'/
290      1          '       ', 6X,'X',11X,'Y',11X,'Z',
291      &                          10X,'X',11X,'Y',11X,'Z')
292   110 format (a,'(',i3,')',6f12.5)
293       
294       enddo
295 cc enddiag
296       do j=1,nbfrag     
297         do i=1,4                                                       
298          bfrag(i,j)=bfrag(i,j)-ishift
299         enddo
300       enddo
301
302       do j=1,nhfrag
303         do i=1,2
304          hfrag(i,j)=hfrag(i,j)-ishift
305         enddo
306       enddo
307
308       return
309       end
310 c---------------------------------------------------------------------------
311       subroutine int_from_cart(lside,lprn)
312       implicit real*8 (a-h,o-z)
313       include 'DIMENSIONS'
314 #ifdef MPI
315       include "mpif.h"
316 #endif 
317       include 'COMMON.LOCAL'
318       include 'COMMON.VAR'
319       include 'COMMON.CHAIN'
320       include 'COMMON.INTERACT'
321       include 'COMMON.IOUNITS'
322       include 'COMMON.GEO'
323       include 'COMMON.NAMES'
324       include 'COMMON.CONTROL'
325       include 'COMMON.SETUP'
326       character*3 seq,atom,res
327       character*80 card
328       dimension sccor(3,20)
329       integer rescode
330       logical lside,lprn
331 #ifdef MPI
332       if(me.eq.king.or..not.out1file)then
333 #endif
334        if (lprn) then 
335         write (iout,'(/a)') 
336      &  'Internal coordinates calculated from crystal structure.'
337         if (lside) then 
338           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
339      & '       Phi','    Dsc_id','       Dsc','     Alpha',
340      & '     Omega'
341         else 
342           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
343      & '       Phi'
344         endif
345        endif
346 #ifdef MPI
347       endif
348 #endif
349       do i=1,nres-1
350         iti=itype(i)
351         if (iti.ne.21 .and. itype(i+1).ne.21 .and. 
352      &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
353           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
354 ctest          stop
355         endif
356         vbld(i+1)=dist(i,i+1)
357         vbld_inv(i+1)=1.0d0/vbld(i+1)
358         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
359         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
360       enddo
361 c      if (unres_pdb) then
362 c        if (itype(1).eq.21) then
363 c          theta(3)=90.0d0*deg2rad
364 c          phi(4)=180.0d0*deg2rad
365 c          vbld(2)=3.8d0
366 c          vbld_inv(2)=1.0d0/vbld(2)
367 c        endif
368 c        if (itype(nres).eq.21) then
369 c          theta(nres)=90.0d0*deg2rad
370 c          phi(nres)=180.0d0*deg2rad
371 c          vbld(nres)=3.8d0
372 c          vbld_inv(nres)=1.0d0/vbld(2)
373 c        endif
374 c      endif
375       if (lside) then
376         do i=2,nres-1
377           do j=1,3
378             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
379      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
380           enddo
381           iti=itype(i)
382           di=dist(i,nres+i)
383           vbld(i+nres)=di
384           if (itype(i).ne.10) then
385             vbld_inv(i+nres)=1.0d0/di
386           else
387             vbld_inv(i+nres)=0.0d0
388           endif
389           if (iti.ne.10) then
390             alph(i)=alpha(nres+i,i,maxres2)
391             omeg(i)=beta(nres+i,i,maxres2,i+1)
392           endif
393           if(me.eq.king.or..not.out1file)then
394            if (lprn)
395      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
396      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
397      &     rad2deg*alph(i),rad2deg*omeg(i)
398           endif
399         enddo
400       else if (lprn) then
401         do i=2,nres
402           iti=itype(i)
403           if(me.eq.king.or..not.out1file)
404      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
405      &     rad2deg*theta(i),rad2deg*phi(i)
406         enddo
407       endif
408       return
409       end
410 c-------------------------------------------------------------------------------
411       subroutine sc_loc_geom(lprn)
412       implicit real*8 (a-h,o-z)
413       include 'DIMENSIONS'
414 #ifdef MPI
415       include "mpif.h"
416 #endif 
417       include 'COMMON.LOCAL'
418       include 'COMMON.VAR'
419       include 'COMMON.CHAIN'
420       include 'COMMON.INTERACT'
421       include 'COMMON.IOUNITS'
422       include 'COMMON.GEO'
423       include 'COMMON.NAMES'
424       include 'COMMON.CONTROL'
425       include 'COMMON.SETUP'
426       double precision x_prime(3),y_prime(3),z_prime(3)
427       logical lprn
428       do i=1,nres-1
429         do j=1,3
430           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
431         enddo
432       enddo
433       do i=2,nres-1
434         if (itype(i).ne.10 .and. itype(i).ne.21) then
435           do j=1,3
436             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
437           enddo
438         else
439           do j=1,3
440             dc_norm(j,i+nres)=0.0d0
441           enddo
442         endif
443       enddo
444       do i=2,nres-1
445         costtab(i+1) =dcos(theta(i+1))
446         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
447         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
448         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
449         cosfac2=0.5d0/(1.0d0+costtab(i+1))
450         cosfac=dsqrt(cosfac2)
451         sinfac2=0.5d0/(1.0d0-costtab(i+1))
452         sinfac=dsqrt(sinfac2)
453         it=itype(i)
454         if (it.ne.10 .and. itype(i).ne.21) then
455 c
456 C  Compute the axes of tghe local cartesian coordinates system; store in
457 c   x_prime, y_prime and z_prime 
458 c
459         do j=1,3
460           x_prime(j) = 0.00
461           y_prime(j) = 0.00
462           z_prime(j) = 0.00
463         enddo
464         do j = 1,3
465           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
466           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
467         enddo
468         call vecpr(x_prime,y_prime,z_prime)
469 c
470 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
471 C to local coordinate system. Store in xx, yy, zz.
472 c
473         xx=0.0d0
474         yy=0.0d0
475         zz=0.0d0
476         do j = 1,3
477           xx = xx + x_prime(j)*dc_norm(j,i+nres)
478           yy = yy + y_prime(j)*dc_norm(j,i+nres)
479           zz = zz + z_prime(j)*dc_norm(j,i+nres)
480         enddo
481
482         xxref(i)=xx
483         yyref(i)=yy
484         zzref(i)=zz
485         else
486         xxref(i)=0.0d0
487         yyref(i)=0.0d0
488         zzref(i)=0.0d0
489         endif
490       enddo
491       if (lprn) then
492         do i=2,nres
493           iti=itype(i)
494 #ifdef MPI
495           if(me.eq.king.or..not.out1file)
496      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
497      &      yyref(i),zzref(i)
498 #else
499           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
500      &     zzref(i)
501 #endif
502         enddo
503       endif
504       return
505       end
506 c---------------------------------------------------------------------------
507       subroutine sccenter(ires,nscat,sccor)
508       implicit real*8 (a-h,o-z)
509       include 'DIMENSIONS'
510       include 'COMMON.CHAIN'
511       dimension sccor(3,20)
512       do j=1,3
513         sccmj=0.0D0
514         do i=1,nscat
515           sccmj=sccmj+sccor(j,i) 
516         enddo
517         dc(j,ires)=sccmj/nscat
518       enddo
519       return
520       end
521 c---------------------------------------------------------------------------
522       subroutine bond_regular
523       implicit real*8 (a-h,o-z)
524       include 'DIMENSIONS'   
525       include 'COMMON.VAR'
526       include 'COMMON.LOCAL'      
527       include 'COMMON.CALC'
528       include 'COMMON.INTERACT'
529       include 'COMMON.CHAIN'
530       do i=1,nres-1
531        vbld(i+1)=vbl
532        vbld_inv(i+1)=1.0d0/vbld(i+1)
533        vbld(i+1+nres)=dsc(itype(i+1))
534        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
535 c       print *,vbld(i+1),vbld(i+1+nres)
536       enddo
537       return
538       end
539