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