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