2 C Read the PDB file and convert the peptide geometry into virtual-chain
4 implicit real*8 (a-h,o-z)
9 include 'COMMON.INTERACT'
10 include 'COMMON.IOUNITS'
12 include 'COMMON.NAMES'
13 include 'COMMON.CONTROL'
14 include 'COMMON.DISTFIT'
15 include 'COMMON.SETUP'
16 character*3 seq,atom,res
19 double precision e1(3),e2(3),e3(3)
27 read (ipdbin,'(a80)',end=10) card
28 if (card(:5).eq.'HELIX') then
31 read(card(22:25),*) hfrag(1,nhfrag)
32 read(card(34:37),*) hfrag(2,nhfrag)
34 if (card(:5).eq.'SHEET') then
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----------------------------------------
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.
54 dc(j,ires+nres)=sccor(j,iii)
57 call sccenter(ires,iii,sccor)
61 read (card(24:26),*) ires
62 read (card(18:20),'(a3)') res
65 if (res.ne.'GLY' .and. res.ne. 'ACE') then
72 if (res.eq.'ACE') then
75 itype(ires)=rescode(ires,res,0)
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)
83 sccor(j,iii)=c(j,ires)
85 else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and.
86 & atom.ne.'N ' .and. atom.ne.'C ') then
88 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
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.
97 dc(j,ires+nres)=sccor(j,iii)
100 call sccenter(ires,iii,sccor)
105 if (itype(nres).ne.10) 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)
117 c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
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)
134 c(j,2*nres)=c(j,nres)
136 if (itype(1).eq.21) 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)
148 c(j,1)=c(j,2)-3.8d0*e2(j)
158 C Calculate internal coordinates.
159 if(me.eq.king.or..not.out1file)then
161 & "Backbone and SC coordinates as read from the PDB"
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)
168 call int_from_cart(.true.,.false.)
169 call sc_loc_geom(.false.)
176 dc(j,i)=c(j,i+1)-c(j,i)
177 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
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)
185 c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
189 C Copy the coordinates to reference coordinates
199 bfrag(i,j)=bfrag(i,j)-ishift
205 hfrag(i,j)=hfrag(i,j)-ishift
211 c---------------------------------------------------------------------------
212 subroutine int_from_cart(lside,lprn)
213 implicit real*8 (a-h,o-z)
216 include 'COMMON.LOCAL'
218 include 'COMMON.CHAIN'
219 include 'COMMON.INTERACT'
220 include 'COMMON.IOUNITS'
222 include 'COMMON.NAMES'
223 include 'COMMON.CONTROL'
224 include 'COMMON.SETUP'
225 character*3 seq,atom,res
227 dimension sccor(3,20)
230 if(me.eq.king.or..not.out1file)then
233 & 'Internal coordinates calculated from crystal structure.'
235 write (iout,'(8a)') ' Res ',' dvb',' Theta',
236 & ' Gamma',' Dsc_id',' Dsc',' Alpha',
239 write (iout,'(4a)') ' Res ',' dvb',' Theta',
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
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) phi(i+1)=beta(i-2,i-1,i,i+1)
255 c if (unres_pdb) then
256 c if (itype(1).eq.21) then
257 c theta(3)=90.0d0*deg2rad
258 c phi(4)=180.0d0*deg2rad
260 c vbld_inv(2)=1.0d0/vbld(2)
262 c if (itype(nres).eq.21) then
263 c theta(nres)=90.0d0*deg2rad
264 c phi(nres)=180.0d0*deg2rad
266 c vbld_inv(nres)=1.0d0/vbld(2)
272 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
273 & +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
277 C 9/29/12 Adam: Correction for zero SC-SC bond length
278 if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
281 if (itype(i).ne.10) then
282 vbld_inv(i+nres)=1.0d0/di
284 vbld_inv(i+nres)=0.0d0
287 alph(i)=alpha(nres+i,i,maxres2)
288 omeg(i)=beta(nres+i,i,maxres2,i+1)
290 if(me.eq.king.or..not.out1file)then
292 & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
293 & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
294 & rad2deg*alph(i),rad2deg*omeg(i)
300 if(me.eq.king.or..not.out1file)
301 & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
302 & rad2deg*theta(i),rad2deg*phi(i)
307 c-------------------------------------------------------------------------------
308 subroutine sc_loc_geom(lprn)
309 implicit real*8 (a-h,o-z)
312 include 'COMMON.LOCAL'
314 include 'COMMON.CHAIN'
315 include 'COMMON.INTERACT'
316 include 'COMMON.IOUNITS'
318 include 'COMMON.NAMES'
319 include 'COMMON.CONTROL'
320 include 'COMMON.SETUP'
321 double precision x_prime(3),y_prime(3),z_prime(3)
325 dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
329 if (itype(i).ne.10) then
331 dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
335 dc_norm(j,i+nres)=0.0d0
340 costtab(i+1) =dcos(theta(i+1))
341 sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
342 cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
343 sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
344 cosfac2=0.5d0/(1.0d0+costtab(i+1))
345 cosfac=dsqrt(cosfac2)
346 sinfac2=0.5d0/(1.0d0-costtab(i+1))
347 sinfac=dsqrt(sinfac2)
351 C Compute the axes of tghe local cartesian coordinates system; store in
352 c x_prime, y_prime and z_prime
360 x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
361 y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
363 call vecpr(x_prime,y_prime,z_prime)
365 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
366 C to local coordinate system. Store in xx, yy, zz.
372 xx = xx + x_prime(j)*dc_norm(j,i+nres)
373 yy = yy + y_prime(j)*dc_norm(j,i+nres)
374 zz = zz + z_prime(j)*dc_norm(j,i+nres)
389 if(me.eq.king.or..not.out1file)
390 & write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
396 c---------------------------------------------------------------------------
397 subroutine sccenter(ires,nscat,sccor)
398 implicit real*8 (a-h,o-z)
400 include 'COMMON.CHAIN'
401 dimension sccor(3,20)
405 sccmj=sccmj+sccor(j,i)
407 dc(j,ires)=sccmj/nscat
411 c---------------------------------------------------------------------------
412 subroutine bond_regular
413 implicit real*8 (a-h,o-z)
416 include 'COMMON.LOCAL'
417 include 'COMMON.CALC'
418 include 'COMMON.INTERACT'
419 include 'COMMON.CHAIN'
422 vbld_inv(i+1)=1.0d0/vbld(i+1)
423 vbld(i+1+nres)=dsc(itype(i+1))
424 vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
425 c print *,vbld(i+1),vbld(i+1+nres)
429 C-----------------------------------------------------------------------------
430 subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
431 implicit real*8 (a-h,o-z)
433 c this subroutine calculates unity vectors of a local reference system
434 c defined by atoms (i2), (i3), and (i4). the x axis is the axis from
435 c atom (i3) to atom (i2), and the xy plane is the plane defined by atoms
436 c (i2), (i3), and (i4). z axis is directed according to the sign of the
437 c vector product (i3)-(i2) and (i3)-(i4). sets fail to .true. if atoms
438 c (i2) and (i3) or (i3) and (i4) coincide or atoms (i2), (i3), and (i4)
439 c form a linear fragment. returns vectors e1, e2, and e3.
441 double precision e1(3),e2(3),e3(3)
442 double precision u(3),z(3)
443 include "COMMON.CHAIN"
444 data coinc /1.0d-13/,align /1.0d-13/
457 if (s1.gt.coinc) goto 2
458 write (iout,1000) i2,i3,i1
461 2 if (s2.gt.coinc) goto 4
462 write(iout,1000) i3,i4,i1
467 v1=z(2)*u(3)-z(3)*u(2)
468 v2=z(3)*u(1)-z(1)*u(3)
469 v3=z(1)*u(2)-z(2)*u(1)
470 anorm=sqrt(v1*v1+v2*v2+v3*v3)
471 if (anorm.gt.align) goto 6
472 write (iout,1010) i2,i3,i4,i1
482 e2(1)=e1(3)*e3(2)-e1(2)*e3(3)
483 e2(2)=e1(1)*e3(3)-e1(3)*e3(1)
484 e2(3)=e1(2)*e3(1)-e1(1)*e3(2)
485 1000 format (/1x,' * * * error - atoms',i4,' and',i4,' coincide.')
486 1010 format (/1x,' * * * error - atoms',2(i4,2h, ),i4,' form a linear')