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