source/cluster/unres/src/chainbuild.f

   1       subroutine chainbuild
   2 C
   3 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
   4 C As of 2/17/95.
   5 C
   6       include 'DIMENSIONS'
   7       include 'COMMON.CHAIN'
   8       include 'COMMON.LOCAL'
   9       include 'COMMON.GEO'
  10       include 'COMMON.VAR'
  11       include 'COMMON.IOUNITS'
  12       include 'COMMON.NAMES'
  13       include 'COMMON.INTERACT'
  14       logical lprn
  15 C Set lprn=.true. for debugging
  16       lprn = .false.
  17 C
  18 C Define the origin and orientation of the coordinate system and locate the
  19 C first three CA's and SC(2).
  20 C
  21       call orig_frame
  22 *
  23 * Build the alpha-carbon chain.
  24 *
  25       do i=4,nres
  26         call locate_next_res(i)
  27       enddo
  28 C
  29 C First and last SC must coincide with the corresponding CA.
  30 C
  31       do j=1,3
  32         dc(j,nres+1)=0.0D0
  33         dc_norm(j,nres+1)=0.0D0
  34         dc(j,nres+nres)=0.0D0
  35         dc_norm(j,nres+nres)=0.0D0
  36         c(j,nres+1)=c(j,1)
  37         c(j,nres+nres)=c(j,nres)
  38       enddo
  39 *
  40 * Temporary diagnosis
  41 *
  42       if (lprn) then
  43
  44       call cartprint
  45       write (iout,'(/a)') 'Recalculated internal coordinates'
  46       do i=2,nres-1
  47         do j=1,3
  48           c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
  49         enddo
  50         be=0.0D0
  51         if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
  52         be1=rad2deg*beta(nres+i,i,maxres2,i+1)
  53         alfai=0.0D0
  54         if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
  55         write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
  56      &  alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
  57       enddo
  58  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
  59
  60       endif
  61
  62       return
  63       end
  64 c-------------------------------------------------------------------------
  65       subroutine orig_frame
  66 C
  67 C Define the origin and orientation of the coordinate system and locate
  68 C the first three atoms.
  69 C
  70       include 'DIMENSIONS'
  71       include 'COMMON.CHAIN'
  72       include 'COMMON.LOCAL'
  73       include 'COMMON.GEO'
  74       include 'COMMON.VAR'
  75       cost=dcos(theta(3))
  76       sint=dsin(theta(3))
  77       t(1,1,1)=-cost
  78       t(1,2,1)=-sint
  79       t(1,3,1)= 0.0D0
  80       t(2,1,1)=-sint
  81       t(2,2,1)= cost
  82       t(2,3,1)= 0.0D0
  83       t(3,1,1)= 0.0D0
  84       t(3,2,1)= 0.0D0
  85       t(3,3,1)= 1.0D0
  86       r(1,1,1)= 1.0D0
  87       r(1,2,1)= 0.0D0
  88       r(1,3,1)= 0.0D0
  89       r(2,1,1)= 0.0D0
  90       r(2,2,1)= 1.0D0
  91       r(2,3,1)= 0.0D0
  92       r(3,1,1)= 0.0D0
  93       r(3,2,1)= 0.0D0
  94       r(3,3,1)= 1.0D0
  95       do i=1,3
  96         do j=1,3
  97           rt(i,j,1)=t(i,j,1)
  98         enddo
  99       enddo
 100       do i=1,3
 101         do j=1,3
 102           prod(i,j,1)=0.0D0
 103           prod(i,j,2)=t(i,j,1)
 104         enddo
 105         prod(i,i,1)=1.0D0
 106       enddo
 107       c(1,1)=0.0D0
 108       c(2,1)=0.0D0
 109       c(3,1)=0.0D0
 110       c(1,2)=vbl
 111       c(2,2)=0.0D0
 112       c(3,2)=0.0D0
 113       dc(1,1)=vbl
 114       dc(2,1)=0.0D0
 115       dc(3,1)=0.0D0
 116       dc_norm(1,1)=1.0D0
 117       dc_norm(2,1)=0.0D0
 118       dc_norm(3,1)=0.0D0
 119       do j=1,3
 120         dc_norm(j,2)=prod(j,1,2)
 121         dc(j,2)=vbl*prod(j,1,2)
 122         c(j,3)=c(j,2)+dc(j,2)
 123       enddo
 124       call locate_side_chain(2)
 125       return
 126       end
 127 c-----------------------------------------------------------------------------
 128       subroutine locate_next_res(i)
 129 C
 130 C Locate CA(i) and SC(i-1)
 131 C
 132       include 'DIMENSIONS'
 133       include 'COMMON.CHAIN'
 134       include 'COMMON.LOCAL'
 135       include 'COMMON.GEO'
 136       include 'COMMON.VAR'
 137       include 'COMMON.IOUNITS'
 138       include 'COMMON.NAMES'
 139       include 'COMMON.INTERACT'
 140 C
 141 C Define the rotation matrices corresponding to CA(i)
 142 C
 143       theti=theta(i)
 144       phii=phi(i)
 145       cost=dcos(theti)
 146       sint=dsin(theti)
 147       cosphi=dcos(phii)
 148       sinphi=dsin(phii)
 149 * Define the matrices of the rotation about the virtual-bond valence angles
 150 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
 151 * program), R(i,j,k), and, the cumulative matrices of rotation RT
 152       t(1,1,i-2)=-cost
 153       t(1,2,i-2)=-sint
 154       t(1,3,i-2)= 0.0D0
 155       t(2,1,i-2)=-sint
 156       t(2,2,i-2)= cost
 157       t(2,3,i-2)= 0.0D0
 158       t(3,1,i-2)= 0.0D0
 159       t(3,2,i-2)= 0.0D0
 160       t(3,3,i-2)= 1.0D0
 161       r(1,1,i-2)= 1.0D0
 162       r(1,2,i-2)= 0.0D0
 163       r(1,3,i-2)= 0.0D0
 164       r(2,1,i-2)= 0.0D0
 165       r(2,2,i-2)=-cosphi
 166       r(2,3,i-2)= sinphi
 167       r(3,1,i-2)= 0.0D0
 168       r(3,2,i-2)= sinphi
 169       r(3,3,i-2)= cosphi
 170       rt(1,1,i-2)=-cost
 171       rt(1,2,i-2)=-sint
 172       rt(1,3,i-2)=0.0D0
 173       rt(2,1,i-2)=sint*cosphi
 174       rt(2,2,i-2)=-cost*cosphi
 175       rt(2,3,i-2)=sinphi
 176       rt(3,1,i-2)=-sint*sinphi
 177       rt(3,2,i-2)=cost*sinphi
 178       rt(3,3,i-2)=cosphi
 179       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
 180       do j=1,3
 181         dc_norm(j,i-1)=prod(j,1,i-1)
 182         dc(j,i-1)=vbl*prod(j,1,i-1)
 183         c(j,i)=c(j,i-1)+dc(j,i-1)
 184       enddo
 185 cd    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
 186 C
 187 C Now calculate the coordinates of SC(i-1)
 188 C
 189       call locate_side_chain(i-1)
 190       return
 191       end
 192 c-----------------------------------------------------------------------------
 193       subroutine locate_side_chain(i)
 194 C
 195 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
 196 C
 197       include 'DIMENSIONS'
 198       include 'COMMON.CHAIN'
 199       include 'COMMON.LOCAL'
 200       include 'COMMON.GEO'
 201       include 'COMMON.VAR'
 202       include 'COMMON.IOUNITS'
 203       include 'COMMON.NAMES'
 204       include 'COMMON.INTERACT'
 205       dimension xx(3)
 206
 207       dsci=dsc(itype(i))
 208       dsci_inv=dsc_inv(itype(i))
 209       alphi=alph(i)
 210       omegi=omeg(i)
 211       cosalphi=dcos(alphi)
 212       sinalphi=dsin(alphi)
 213       cosomegi=dcos(omegi)
 214       sinomegi=dsin(omegi)
 215       xp= dsci*cosalphi
 216       yp= dsci*sinalphi*cosomegi
 217       zp=-dsci*sinalphi*sinomegi
 218 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
 219 * X-axis aligned with the vector DC(*,i)
 220       theta2=pi-0.5D0*theta(i+1)
 221       cost2=dcos(theta2)
 222       sint2=dsin(theta2)
 223       xx(1)= xp*cost2+yp*sint2
 224       xx(2)=-xp*sint2+yp*cost2
 225       xx(3)= zp
 226 cd    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
 227 cd   &   xp,yp,zp,(xx(k),k=1,3)
 228       do j=1,3
 229         xloc(j,i)=xx(j)
 230       enddo
 231 * Bring the SC vectors to the common coordinate system.
 232       xx(1)=xloc(1,i)
 233       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
 234       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
 235       do j=1,3
 236         xrot(j,i)=xx(j)
 237       enddo
 238       do j=1,3
 239         rj=0.0D0
 240         do k=1,3
 241           rj=rj+prod(j,k,i-1)*xx(k)
 242         enddo
 243         dc(j,nres+i)=rj
 244         dc_norm(j,nres+i)=rj*dsci_inv
 245         c(j,nres+i)=c(j,i)+rj
 246       enddo
 247       return
 248       end