source/wham/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 '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 #ifdef OSF
 150       theti=theta(i)
 151       icrc=0
 152       call proc_proc(theti,icrc)
 153       if(icrc.eq.1)theti=100.0
 154       phii=phi(i)
 155       icrc=0
 156       call proc_proc(phii,icrc)
 157       if(icrc.eq.1)phii=180.0
 158 #else
 159       theti=theta(i)
 160       phii=phi(i)
 161 #endif
 162       cost=dcos(theti)
 163       sint=dsin(theti)
 164       cosphi=dcos(phii)
 165       sinphi=dsin(phii)
 166 * Define the matrices of the rotation about the virtual-bond valence angles
 167 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
 168 * program), R(i,j,k), and, the cumulative matrices of rotation RT
 169       t(1,1,i-2)=-cost
 170       t(1,2,i-2)=-sint
 171       t(1,3,i-2)= 0.0D0
 172       t(2,1,i-2)=-sint
 173       t(2,2,i-2)= cost
 174       t(2,3,i-2)= 0.0D0
 175       t(3,1,i-2)= 0.0D0
 176       t(3,2,i-2)= 0.0D0
 177       t(3,3,i-2)= 1.0D0
 178       r(1,1,i-2)= 1.0D0
 179       r(1,2,i-2)= 0.0D0
 180       r(1,3,i-2)= 0.0D0
 181       r(2,1,i-2)= 0.0D0
 182       r(2,2,i-2)=-cosphi
 183       r(2,3,i-2)= sinphi
 184       r(3,1,i-2)= 0.0D0
 185       r(3,2,i-2)= sinphi
 186       r(3,3,i-2)= cosphi
 187       rt(1,1,i-2)=-cost
 188       rt(1,2,i-2)=-sint
 189       rt(1,3,i-2)=0.0D0
 190       rt(2,1,i-2)=sint*cosphi
 191       rt(2,2,i-2)=-cost*cosphi
 192       rt(2,3,i-2)=sinphi
 193       rt(3,1,i-2)=-sint*sinphi
 194       rt(3,2,i-2)=cost*sinphi
 195       rt(3,3,i-2)=cosphi
 196       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
 197       do j=1,3
 198         dc_norm(j,i-1)=prod(j,1,i-1)
 199         dc(j,i-1)=vbld(i)*prod(j,1,i-1)
 200         c(j,i)=c(j,i-1)+dc(j,i-1)
 201       enddo
 202 cd    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
 203 C
 204 C Now calculate the coordinates of SC(i-1)
 205 C
 206       call locate_side_chain(i-1)
 207       return
 208       end
 209 c-----------------------------------------------------------------------------
 210       subroutine locate_side_chain(i)
 211 C
 212 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
 213 C
 214       implicit real*8 (a-h,o-z)
 215       include 'DIMENSIONS'
 216       include 'DIMENSIONS.ZSCOPT'
 217       include 'COMMON.CHAIN'
 218       include 'COMMON.LOCAL'
 219       include 'COMMON.GEO'
 220       include 'COMMON.VAR'
 221       include 'COMMON.IOUNITS'
 222       include 'COMMON.NAMES'
 223       include 'COMMON.INTERACT'
 224       dimension xx(3)
 225
 226 c      dsci=dsc(itype(i))
 227 c      dsci_inv=dsc_inv(itype(i))
 228       dsci=vbld(i+nres)
 229       dsci_inv=vbld_inv(i+nres)
 230 #ifdef OSF
 231       alphi=alph(i)
 232       omegi=omeg(i)
 233 c detecting NaNQ
 234       icrc=0
 235       call proc_proc(alphi,icrc)
 236       if(icrc.eq.1)alphi=100.0
 237       icrc=0
 238       call proc_proc(omegi,icrc)
 239       if(icrc.eq.1)omegi=-100.0
 240 #else
 241       alphi=alph(i)
 242       omegi=omeg(i)
 243 #endif
 244       cosalphi=dcos(alphi)
 245       sinalphi=dsin(alphi)
 246       cosomegi=dcos(omegi)
 247       sinomegi=dsin(omegi)
 248       xp= dsci*cosalphi
 249       yp= dsci*sinalphi*cosomegi
 250       zp=-dsci*sinalphi*sinomegi
 251 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
 252 * X-axis aligned with the vector DC(*,i)
 253       theta2=pi-0.5D0*theta(i+1)
 254       cost2=dcos(theta2)
 255       sint2=dsin(theta2)
 256       xx(1)= xp*cost2+yp*sint2
 257       xx(2)=-xp*sint2+yp*cost2
 258       xx(3)= zp
 259 cd    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
 260 cd   &   xp,yp,zp,(xx(k),k=1,3)
 261       do j=1,3
 262         xloc(j,i)=xx(j)
 263       enddo
 264 * Bring the SC vectors to the common coordinate system.
 265       xx(1)=xloc(1,i)
 266       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
 267       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
 268       do j=1,3
 269         xrot(j,i)=xx(j)
 270       enddo
 271       do j=1,3
 272         rj=0.0D0
 273         do k=1,3
 274           rj=rj+prod(j,k,i-1)*xx(k)
 275         enddo
 276         dc(j,nres+i)=rj
 277         dc_norm(j,nres+i)=rj*dsci_inv
 278         c(j,nres+i)=c(j,i)+rj
 279       enddo
 280       return
 281       end