source/unres/src_MD-M/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       double precision e1(3),e2(3),e3(3)
  16       logical lprn,perbox,fail
  17 C Set lprn=.true. for debugging
  18       lprn = .false.
  19       perbox=.false.
  20       fail=.false.
  21       print *, 'enter chainbuild'
  22       call chainbuild_cart
  23       return
  24       end
  25 #ifdef DEBUG
  26       if (perbox) then
  27       cost=dcos(theta(3))
  28       sint=dsin(theta(3))
  29       print *,'before refsys'
  30       call refsys(2,3,4,e1,e2,e3,fail)
  31       print *,'after refsys'
  32           if (fail) then
  33             e2(1)=0.0d0
  34             e2(2)=1.0d0
  35             e2(3)=0.0d0
  36           endif
  37       dc(1,0)=c(1,1)
  38       dc(2,0)=c(2,1)
  39       dc(3,0)=c(3,1)
  40       print *,'dc',dc(1,0),dc(2,0),dc(3,0)
  41       dc(1,1)=c(1,2)-c(1,1)
  42       dc(2,1)=c(2,2)-c(2,1)
  43       dc(3,1)=c(3,2)-c(3,1)
  44       dc(1,2)=c(1,3)-c(1,2)
  45       dc(2,2)=c(2,3)-c(2,2)
  46       dc(3,2)=c(3,3)-c(3,2)
  47       t(1,1,1)=e1(1)
  48       t(1,2,1)=e1(2)
  49       t(1,3,1)=e1(3)
  50       t(2,1,1)=e2(1)
  51       t(2,2,1)=e2(2)
  52       t(2,3,1)=e2(3)
  53       t(3,1,1)=e3(1)
  54       t(3,2,1)=e3(2)
  55       t(3,3,1)=e3(3)
  56       veclen=0.0d0
  57       do i=1,3
  58        veclen=veclen+(c(i,2)-c(i,1))**2
  59       enddo
  60       veclen=sqrt(veclen)
  61       r(1,1,1)= 1.0D0
  62       r(1,2,1)= 0.0D0
  63       r(1,3,1)= 0.0D0
  64       r(2,1,1)= 0.0D0
  65       r(2,2,1)= 1.0D0
  66       r(2,3,1)= 0.0D0
  67       r(3,1,1)= 0.0D0
  68       r(3,2,1)= 0.0D0
  69       r(3,3,1)= 1.0D0
  70       do i=1,3
  71         do j=1,3
  72           rt(i,j,1)=t(i,j,1)
  73         enddo
  74       enddo
  75       do i=1,3
  76         do j=1,3
  77           prod(i,j,1)=0.0D0
  78           prod(i,j,2)=t(i,j,1)
  79         enddo
  80         prod(i,i,1)=1.0D0
  81       enddo
  82         call locate_side_chain(2)
  83       do i=4,nres
  84 #ifdef OSF
  85       theti=theta(i)
  86       if (theti.ne.theti) theti=100.0
  87       phii=phi(i)
  88       if (phii.ne.phii) phii=180.0
  89 #else
  90       theti=theta(i)
  91       phii=phi(i)
  92 #endif
  93       cost=dcos(theti)
  94       sint=dsin(theti)
  95       cosphi=dcos(phii)
  96       sinphi=dsin(phii)
  97 * Define the matrices of the rotation about the virtual-bond valence angles
  98 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
  99 * program), R(i,j,k), and, the cumulative matrices of rotation RT
 100       t(1,1,i-2)=-cost
 101       t(1,2,i-2)=-sint
 102       t(1,3,i-2)= 0.0D0
 103       t(2,1,i-2)=-sint
 104       t(2,2,i-2)= cost
 105       t(2,3,i-2)= 0.0D0
 106       t(3,1,i-2)= 0.0D0
 107       t(3,2,i-2)= 0.0D0
 108       t(3,3,i-2)= 1.0D0
 109       r(1,1,i-2)= 1.0D0
 110       r(1,2,i-2)= 0.0D0
 111       r(1,3,i-2)= 0.0D0
 112       r(2,1,i-2)= 0.0D0
 113       r(2,2,i-2)=-cosphi
 114       r(2,3,i-2)= sinphi
 115       r(3,1,i-2)= 0.0D0
 116       r(3,2,i-2)= sinphi
 117       r(3,3,i-2)= cosphi
 118       rt(1,1,i-2)=-cost
 119       rt(1,2,i-2)=-sint
 120       rt(1,3,i-2)=0.0D0
 121       rt(2,1,i-2)=sint*cosphi
 122       rt(2,2,i-2)=-cost*cosphi
 123       rt(2,3,i-2)=sinphi
 124       rt(3,1,i-2)=-sint*sinphi
 125       rt(3,2,i-2)=cost*sinphi
 126       rt(3,3,i-2)=cosphi
 127         call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
 128       do j=1,3
 129         dc_norm(j,i-1)=prod(j,1,i-1)
 130         dc(j,i-1)=vbld(i)*prod(j,1,i-1)
 131       enddo
 132         call locate_side_chain(i-1)
 133        enddo
 134       else
 135 C
 136 C Define the origin and orientation of the coordinate system and locate the
 137 C first three CA's and SC(2).
 138 C
 139       call orig_frame
 140 *
 141 * Build the alpha-carbon chain.
 142 *
 143       do i=4,nres
 144         call locate_next_res(i)
 145       enddo
 146 C
 147 C First and last SC must coincide with the corresponding CA.
 148 C
 149       do j=1,3
 150         dc(j,nres+1)=0.0D0
 151         dc_norm(j,nres+1)=0.0D0
 152         dc(j,nres+nres)=0.0D0
 153         dc_norm(j,nres+nres)=0.0D0
 154         c(j,nres+1)=c(j,1)
 155         c(j,nres+nres)=c(j,nres)
 156       enddo
 157 *
 158 * Temporary diagnosis
 159 *
 160       if (lprn) then
 161
 162       call cartprint
 163       write (iout,'(/a)') 'Recalculated internal coordinates'
 164       do i=2,nres-1
 165         do j=1,3
 166           c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 167         enddo
 168         be=0.0D0
 169         if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
 170         be1=rad2deg*beta(nres+i,i,maxres2,i+1)
 171         alfai=0.0D0
 172         if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
 173         write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
 174      &  alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
 175       enddo
 176  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
 177
 178       endif
 179       endif
 180       return
 181       end
 182 #endif
 183 c-------------------------------------------------------------------------
 184       subroutine orig_frame
 185 C
 186 C Define the origin and orientation of the coordinate system and locate
 187 C the first three atoms.
 188 C
 189       implicit real*8 (a-h,o-z)
 190       include 'DIMENSIONS'
 191       include 'COMMON.CHAIN'
 192       include 'COMMON.LOCAL'
 193       include 'COMMON.GEO'
 194       include 'COMMON.VAR'
 195       cost=dcos(theta(3))
 196       sint=dsin(theta(3))
 197       t(1,1,1)=-cost
 198       t(1,2,1)=-sint
 199       t(1,3,1)= 0.0D0
 200       t(2,1,1)=-sint
 201       t(2,2,1)= cost
 202       t(2,3,1)= 0.0D0
 203       t(3,1,1)= 0.0D0
 204       t(3,2,1)= 0.0D0
 205       t(3,3,1)= 1.0D0
 206       r(1,1,1)= 1.0D0
 207       r(1,2,1)= 0.0D0
 208       r(1,3,1)= 0.0D0
 209       r(2,1,1)= 0.0D0
 210       r(2,2,1)= 1.0D0
 211       r(2,3,1)= 0.0D0
 212       r(3,1,1)= 0.0D0
 213       r(3,2,1)= 0.0D0
 214       r(3,3,1)= 1.0D0
 215       do i=1,3
 216         do j=1,3
 217           rt(i,j,1)=t(i,j,1)
 218         enddo
 219       enddo
 220       do i=1,3
 221         do j=1,3
 222           prod(i,j,1)=0.0D0
 223           prod(i,j,2)=t(i,j,1)
 224         enddo
 225         prod(i,i,1)=1.0D0
 226       enddo
 227       c(1,1)=0.0D0
 228       c(2,1)=0.0D0
 229       c(3,1)=0.0D0
 230       c(1,2)=vbld(2)
 231       c(2,2)=0.0D0
 232       c(3,2)=0.0D0
 233       dc(1,0)=0.0d0
 234       dc(2,0)=0.0D0
 235       dc(3,0)=0.0D0
 236       dc(1,1)=vbld(2)
 237       dc(2,1)=0.0D0
 238       dc(3,1)=0.0D0
 239       dc_norm(1,0)=0.0D0
 240       dc_norm(2,0)=0.0D0
 241       dc_norm(3,0)=0.0D0
 242       dc_norm(1,1)=1.0D0
 243       dc_norm(2,1)=0.0D0
 244       dc_norm(3,1)=0.0D0
 245       do j=1,3
 246         dc_norm(j,2)=prod(j,1,2)
 247         dc(j,2)=vbld(3)*prod(j,1,2)
 248         c(j,3)=c(j,2)+dc(j,2)
 249       enddo
 250       call locate_side_chain(2)
 251       return
 252       end
 253 c-----------------------------------------------------------------------------
 254       subroutine locate_next_res(i)
 255 C
 256 C Locate CA(i) and SC(i-1)
 257 C
 258       implicit real*8 (a-h,o-z)
 259       include 'DIMENSIONS'
 260       include 'COMMON.CHAIN'
 261       include 'COMMON.LOCAL'
 262       include 'COMMON.GEO'
 263       include 'COMMON.VAR'
 264       include 'COMMON.IOUNITS'
 265       include 'COMMON.NAMES'
 266       include 'COMMON.INTERACT'
 267 C
 268 C Define the rotation matrices corresponding to CA(i)
 269 C
 270 #ifdef OSF
 271       theti=theta(i)
 272       if (theti.ne.theti) theti=100.0
 273       phii=phi(i)
 274       if (phii.ne.phii) phii=180.0
 275 #else
 276       theti=theta(i)
 277       phii=phi(i)
 278 #endif
 279       cost=dcos(theti)
 280       sint=dsin(theti)
 281       cosphi=dcos(phii)
 282       sinphi=dsin(phii)
 283 * Define the matrices of the rotation about the virtual-bond valence angles
 284 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
 285 * program), R(i,j,k), and, the cumulative matrices of rotation RT
 286       t(1,1,i-2)=-cost
 287       t(1,2,i-2)=-sint
 288       t(1,3,i-2)= 0.0D0
 289       t(2,1,i-2)=-sint
 290       t(2,2,i-2)= cost
 291       t(2,3,i-2)= 0.0D0
 292       t(3,1,i-2)= 0.0D0
 293       t(3,2,i-2)= 0.0D0
 294       t(3,3,i-2)= 1.0D0
 295       r(1,1,i-2)= 1.0D0
 296       r(1,2,i-2)= 0.0D0
 297       r(1,3,i-2)= 0.0D0
 298       r(2,1,i-2)= 0.0D0
 299       r(2,2,i-2)=-cosphi
 300       r(2,3,i-2)= sinphi
 301       r(3,1,i-2)= 0.0D0
 302       r(3,2,i-2)= sinphi
 303       r(3,3,i-2)= cosphi
 304       rt(1,1,i-2)=-cost
 305       rt(1,2,i-2)=-sint
 306       rt(1,3,i-2)=0.0D0
 307       rt(2,1,i-2)=sint*cosphi
 308       rt(2,2,i-2)=-cost*cosphi
 309       rt(2,3,i-2)=sinphi
 310       rt(3,1,i-2)=-sint*sinphi
 311       rt(3,2,i-2)=cost*sinphi
 312       rt(3,3,i-2)=cosphi
 313       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
 314       do j=1,3
 315         dc_norm(j,i-1)=prod(j,1,i-1)
 316         dc(j,i-1)=vbld(i)*prod(j,1,i-1)
 317         c(j,i)=c(j,i-1)+dc(j,i-1)
 318       enddo
 319 cd    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
 320 C
 321 C Now calculate the coordinates of SC(i-1)
 322 C
 323       call locate_side_chain(i-1)
 324       return
 325       end
 326 c-----------------------------------------------------------------------------
 327       subroutine locate_side_chain(i)
 328 C
 329 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
 330 C
 331       implicit real*8 (a-h,o-z)
 332       include 'DIMENSIONS'
 333       include 'COMMON.CHAIN'
 334       include 'COMMON.LOCAL'
 335       include 'COMMON.GEO'
 336       include 'COMMON.VAR'
 337       include 'COMMON.IOUNITS'
 338       include 'COMMON.NAMES'
 339       include 'COMMON.INTERACT'
 340       dimension xx(3)
 341
 342 c      dsci=dsc(itype(i))
 343 c      dsci_inv=dsc_inv(itype(i))
 344       dsci=vbld(i+nres)
 345       dsci_inv=vbld_inv(i+nres)
 346 #ifdef OSF
 347       alphi=alph(i)
 348       omegi=omeg(i)
 349       if (alphi.ne.alphi) alphi=100.0
 350       if (omegi.ne.omegi) omegi=-100.0
 351 #else
 352       alphi=alph(i)
 353       omegi=omeg(i)
 354 #endif
 355       cosalphi=dcos(alphi)
 356       sinalphi=dsin(alphi)
 357       cosomegi=dcos(omegi)
 358       sinomegi=dsin(omegi)
 359       xp= dsci*cosalphi
 360       yp= dsci*sinalphi*cosomegi
 361       zp=-dsci*sinalphi*sinomegi
 362 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
 363 * X-axis aligned with the vector DC(*,i)
 364       theta2=pi-0.5D0*theta(i+1)
 365       cost2=dcos(theta2)
 366       sint2=dsin(theta2)
 367       xx(1)= xp*cost2+yp*sint2
 368       xx(2)=-xp*sint2+yp*cost2
 369       xx(3)= zp
 370 cd    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
 371 cd   &   xp,yp,zp,(xx(k),k=1,3)
 372       do j=1,3
 373         xloc(j,i)=xx(j)
 374       enddo
 375 * Bring the SC vectors to the common coordinate system.
 376       xx(1)=xloc(1,i)
 377       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
 378       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
 379       do j=1,3
 380         xrot(j,i)=xx(j)
 381       enddo
 382       do j=1,3
 383         rj=0.0D0
 384         do k=1,3
 385           rj=rj+prod(j,k,i-1)*xx(k)
 386         enddo
 387         dc(j,nres+i)=rj
 388         dc_norm(j,nres+i)=rj*dsci_inv
 389         c(j,nres+i)=c(j,i)+rj
 390       enddo
 391       return
 392       end