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