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