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