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