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