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       call chainbuild_cart
  22       return
  23       end
  24 C#ifdef DEBUG
  25 C      if (perbox) then
  26 C first three CA's and SC(2).
  27 C
  28       subroutine chainbuild_extconf
  29 C
  30 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
  31 C As of 2/17/95.
  32 C
  33       implicit real*8 (a-h,o-z)
  34       include 'DIMENSIONS'
  35       include 'COMMON.CHAIN'
  36       include 'COMMON.LOCAL'
  37       include 'COMMON.GEO'
  38       include 'COMMON.VAR'
  39       include 'COMMON.IOUNITS'
  40       include 'COMMON.NAMES'
  41       include 'COMMON.INTERACT'
  42       double precision e1(3),e2(3),e3(3)
  43       logical lprn,perbox,fail
  44
  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 C      endif
  85       endif
  86       return
  87       end
  88 C#endif
  89 c-------------------------------------------------------------------------
  90       subroutine orig_frame
  91 C
  92 C Define the origin and orientation of the coordinate system and locate
  93 C the first three atoms.
  94 C
  95       implicit real*8 (a-h,o-z)
  96       include 'DIMENSIONS'
  97       include 'COMMON.CHAIN'
  98       include 'COMMON.LOCAL'
  99       include 'COMMON.GEO'
 100       include 'COMMON.VAR'
 101       cost=dcos(theta(3))
 102       sint=dsin(theta(3))
 103       t(1,1,1)=-cost
 104       t(1,2,1)=-sint
 105       t(1,3,1)= 0.0D0
 106       t(2,1,1)=-sint
 107       t(2,2,1)= cost
 108       t(2,3,1)= 0.0D0
 109       t(3,1,1)= 0.0D0
 110       t(3,2,1)= 0.0D0
 111       t(3,3,1)= 1.0D0
 112       r(1,1,1)= 1.0D0
 113       r(1,2,1)= 0.0D0
 114       r(1,3,1)= 0.0D0
 115       r(2,1,1)= 0.0D0
 116       r(2,2,1)= 1.0D0
 117       r(2,3,1)= 0.0D0
 118       r(3,1,1)= 0.0D0
 119       r(3,2,1)= 0.0D0
 120       r(3,3,1)= 1.0D0
 121       do i=1,3
 122         do j=1,3
 123           rt(i,j,1)=t(i,j,1)
 124         enddo
 125       enddo
 126       do i=1,3
 127         do j=1,3
 128           prod(i,j,1)=0.0D0
 129           prod(i,j,2)=t(i,j,1)
 130         enddo
 131         prod(i,i,1)=1.0D0
 132       enddo
 133       c(1,1)=0.0D0
 134       c(2,1)=0.0D0
 135       c(3,1)=0.0D0
 136       c(1,2)=vbld(2)
 137       c(2,2)=0.0D0
 138       c(3,2)=0.0D0
 139       dc(1,0)=0.0d0
 140       dc(2,0)=0.0D0
 141       dc(3,0)=0.0D0
 142       dc(1,1)=vbld(2)
 143       dc(2,1)=0.0D0
 144       dc(3,1)=0.0D0
 145       dc_norm(1,0)=0.0D0
 146       dc_norm(2,0)=0.0D0
 147       dc_norm(3,0)=0.0D0
 148       dc_norm(1,1)=1.0D0
 149       dc_norm(2,1)=0.0D0
 150       dc_norm(3,1)=0.0D0
 151       do j=1,3
 152         dc_norm(j,2)=prod(j,1,2)
 153         dc(j,2)=vbld(3)*prod(j,1,2)
 154         c(j,3)=c(j,2)+dc(j,2)
 155       enddo
 156       call locate_side_chain(2)
 157       return
 158       end
 159 c-----------------------------------------------------------------------------
 160       subroutine locate_next_res(i)
 161 C
 162 C Locate CA(i) and SC(i-1)
 163 C
 164       implicit real*8 (a-h,o-z)
 165       include 'DIMENSIONS'
 166       include 'COMMON.CHAIN'
 167       include 'COMMON.LOCAL'
 168       include 'COMMON.GEO'
 169       include 'COMMON.VAR'
 170       include 'COMMON.IOUNITS'
 171       include 'COMMON.NAMES'
 172       include 'COMMON.INTERACT'
 173 C
 174 C Define the rotation matrices corresponding to CA(i)
 175 C
 176 #ifdef OSF
 177       theti=theta(i)
 178       if (theti.ne.theti) theti=100.0
 179       phii=phi(i)
 180       if (phii.ne.phii) phii=180.0
 181 #else
 182       theti=theta(i)
 183       phii=phi(i)
 184 #endif
 185       cost=dcos(theti)
 186       sint=dsin(theti)
 187       cosphi=dcos(phii)
 188       sinphi=dsin(phii)
 189 * Define the matrices of the rotation about the virtual-bond valence angles
 190 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
 191 * program), R(i,j,k), and, the cumulative matrices of rotation RT
 192       t(1,1,i-2)=-cost
 193       t(1,2,i-2)=-sint
 194       t(1,3,i-2)= 0.0D0
 195       t(2,1,i-2)=-sint
 196       t(2,2,i-2)= cost
 197       t(2,3,i-2)= 0.0D0
 198       t(3,1,i-2)= 0.0D0
 199       t(3,2,i-2)= 0.0D0
 200       t(3,3,i-2)= 1.0D0
 201       r(1,1,i-2)= 1.0D0
 202       r(1,2,i-2)= 0.0D0
 203       r(1,3,i-2)= 0.0D0
 204       r(2,1,i-2)= 0.0D0
 205       r(2,2,i-2)=-cosphi
 206       r(2,3,i-2)= sinphi
 207       r(3,1,i-2)= 0.0D0
 208       r(3,2,i-2)= sinphi
 209       r(3,3,i-2)= cosphi
 210       rt(1,1,i-2)=-cost
 211       rt(1,2,i-2)=-sint
 212       rt(1,3,i-2)=0.0D0
 213       rt(2,1,i-2)=sint*cosphi
 214       rt(2,2,i-2)=-cost*cosphi
 215       rt(2,3,i-2)=sinphi
 216       rt(3,1,i-2)=-sint*sinphi
 217       rt(3,2,i-2)=cost*sinphi
 218       rt(3,3,i-2)=cosphi
 219       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
 220       do j=1,3
 221         dc_norm(j,i-1)=prod(j,1,i-1)
 222         dc(j,i-1)=vbld(i)*prod(j,1,i-1)
 223         c(j,i)=c(j,i-1)+dc(j,i-1)
 224       enddo
 225 cd    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
 226 C
 227 C Now calculate the coordinates of SC(i-1)
 228 C
 229       call locate_side_chain(i-1)
 230       return
 231       end
 232 c-----------------------------------------------------------------------------
 233       subroutine locate_side_chain(i)
 234 C
 235 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
 236 C
 237       implicit real*8 (a-h,o-z)
 238       include 'DIMENSIONS'
 239       include 'COMMON.CHAIN'
 240       include 'COMMON.LOCAL'
 241       include 'COMMON.GEO'
 242       include 'COMMON.VAR'
 243       include 'COMMON.IOUNITS'
 244       include 'COMMON.NAMES'
 245       include 'COMMON.INTERACT'
 246       dimension xx(3)
 247
 248 c      dsci=dsc(itype(i))
 249 c      dsci_inv=dsc_inv(itype(i))
 250       dsci=vbld(i+nres)
 251       dsci_inv=vbld_inv(i+nres)
 252 #ifdef OSF
 253       alphi=alph(i)
 254       omegi=omeg(i)
 255       if (alphi.ne.alphi) alphi=100.0
 256       if (omegi.ne.omegi) omegi=-100.0
 257 #else
 258       alphi=alph(i)
 259       omegi=omeg(i)
 260 #endif
 261       cosalphi=dcos(alphi)
 262       sinalphi=dsin(alphi)
 263       cosomegi=dcos(omegi)
 264       sinomegi=dsin(omegi)
 265       xp= dsci*cosalphi
 266       yp= dsci*sinalphi*cosomegi
 267       zp=-dsci*sinalphi*sinomegi
 268 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
 269 * X-axis aligned with the vector DC(*,i)
 270       theta2=pi-0.5D0*theta(i+1)
 271       cost2=dcos(theta2)
 272       sint2=dsin(theta2)
 273       xx(1)= xp*cost2+yp*sint2
 274       xx(2)=-xp*sint2+yp*cost2
 275       xx(3)= zp
 276 cd    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
 277 cd   &   xp,yp,zp,(xx(k),k=1,3)
 278       do j=1,3
 279         xloc(j,i)=xx(j)
 280       enddo
 281 * Bring the SC vectors to the common coordinate system.
 282       xx(1)=xloc(1,i)
 283       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
 284       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
 285       do j=1,3
 286         xrot(j,i)=xx(j)
 287       enddo
 288       do j=1,3
 289         rj=0.0D0
 290         do k=1,3
 291           rj=rj+prod(j,k,i-1)*xx(k)
 292         enddo
 293         dc(j,nres+i)=rj
 294         dc_norm(j,nres+i)=rj*dsci_inv
 295         c(j,nres+i)=c(j,i)+rj
 296       enddo
 297       return
 298       end
 299 c------------------------------------------
 300       subroutine returnbox
 301       include 'DIMENSIONS'
 302 #ifdef MPI
 303       include 'mpif.h'
 304 #endif
 305       include 'COMMON.CONTROL'
 306       include 'COMMON.VAR'
 307       include 'COMMON.MD'
 308 #ifndef LANG0
 309       include 'COMMON.LANGEVIN'
 310 #else
 311       include 'COMMON.LANGEVIN.lang0'
 312 #endif
 313       include 'COMMON.CHAIN'
 314       include 'COMMON.DERIV'
 315       include 'COMMON.GEO'
 316       include 'COMMON.LOCAL'
 317       include 'COMMON.INTERACT'
 318       include 'COMMON.IOUNITS'
 319       include 'COMMON.NAMES'
 320       include 'COMMON.TIME1'
 321       include 'COMMON.REMD'
 322       include 'COMMON.SETUP'
 323       include 'COMMON.MUCA'
 324       include 'COMMON.HAIRPIN'
 325 C change suggested by Ana - begin
 326       integer allareout
 327 C change suggested by Ana - end
 328         j=1
 329         chain_beg=1
 330 C        do i=1,nres
 331 C       write(*,*) 'initial', i,j,c(j,i)
 332 C        enddo
 333 C change suggested by Ana - begin
 334         allareout=1
 335 C change suggested by Ana -end
 336         do i=1,nres-1
 337          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 338           chain_end=i
 339           if (allareout.eq.1) then
 340             ireturnval=int(c(j,i)/boxxsize)
 341             if (c(j,i).le.0) ireturnval=ireturnval-1
 342             do k=chain_beg,chain_end
 343               c(j,k)=c(j,k)-ireturnval*boxxsize
 344               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
 345             enddo
 346 C Suggested by Ana
 347             if (chain_beg.eq.1)
 348      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 349 C Suggested by Ana -end
 350            endif
 351            chain_beg=i+1
 352            allareout=1
 353          else
 354           if (int(c(j,i)/boxxsize).eq.0) allareout=0
 355          endif
 356         enddo
 357          if (allareout.eq.1) then
 358             ireturnval=int(c(j,i)/boxxsize)
 359             if (c(j,i).le.0) ireturnval=ireturnval-1
 360             do k=chain_beg,nres
 361               c(j,k)=c(j,k)-ireturnval*boxxsize
 362               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
 363             enddo
 364           endif
 365 C NO JUMP
 366 C        do i=1,nres
 367 C        write(*,*) 'befor no jump', i,j,c(j,i)
 368 C        enddo
 369         nojumpval=0
 370         do i=2,nres
 371            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 372              difference=abs(c(j,i-1)-c(j,i))
 373 C             print *,'diff', difference
 374              if (difference.gt.boxxsize/2.0) then
 375                 if (c(j,i-1).gt.c(j,i)) then
 376                   nojumpval=1
 377                  else
 378                    nojumpval=-1
 379                  endif
 380               else
 381               nojumpval=0
 382               endif
 383               endif
 384               c(j,i)=c(j,i)+nojumpval*boxxsize
 385               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
 386          enddo
 387        nojumpval=0
 388         do i=2,nres
 389            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 390              difference=abs(c(j,i-1)-c(j,i))
 391              if (difference.gt.boxxsize/2.0) then
 392                 if (c(j,i-1).gt.c(j,i)) then
 393                   nojumpval=1
 394                  else
 395                    nojumpval=-1
 396                  endif
 397               else
 398               nojumpval=0
 399               endif
 400              endif
 401               c(j,i)=c(j,i)+nojumpval*boxxsize
 402               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
 403          enddo
 404
 405 C        do i=1,nres
 406 C        write(*,*) 'after no jump', i,j,c(j,i)
 407 C        enddo
 408
 409 C NOW Y dimension
 410 C suggesed by Ana begins
 411         allareout=1
 412 C suggested by Ana ends
 413         j=2
 414         chain_beg=1
 415         do i=1,nres-1
 416          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 417           chain_end=i
 418           if (allareout.eq.1) then
 419             ireturnval=int(c(j,i)/boxysize)
 420             if (c(j,i).le.0) ireturnval=ireturnval-1
 421             do k=chain_beg,chain_end
 422               c(j,k)=c(j,k)-ireturnval*boxysize
 423              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
 424             enddo
 425 C Suggested by Ana
 426             if (chain_beg.eq.1)
 427      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 428 C Suggested by Ana -end
 429            endif
 430            chain_beg=i+1
 431            allareout=1
 432          else
 433           if (int(c(j,i)/boxysize).eq.0) allareout=0
 434          endif
 435         enddo
 436          if (allareout.eq.1) then
 437             ireturnval=int(c(j,i)/boxysize)
 438             if (c(j,i).le.0) ireturnval=ireturnval-1
 439             do k=chain_beg,nres
 440               c(j,k)=c(j,k)-ireturnval*boxysize
 441               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
 442             enddo
 443           endif
 444         nojumpval=0
 445         do i=2,nres
 446            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 447              difference=abs(c(j,i-1)-c(j,i))
 448              if (difference.gt.boxysize/2.0) then
 449                 if (c(j,i-1).gt.c(j,i)) then
 450                   nojumpval=1
 451                  else
 452                    nojumpval=-1
 453                  endif
 454              else
 455               nojumpval=0
 456               endif
 457            endif
 458               c(j,i)=c(j,i)+nojumpval*boxysize
 459               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
 460          enddo
 461       nojumpval=0
 462         do i=2,nres
 463            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 464              difference=abs(c(j,i-1)-c(j,i))
 465              if (difference.gt.boxysize/2.0) then
 466                 if (c(j,i-1).gt.c(j,i)) then
 467                   nojumpval=1
 468                  else
 469                    nojumpval=-1
 470                  endif
 471               else
 472               nojumpval=0
 473               endif
 474             endif
 475               c(j,i)=c(j,i)+nojumpval*boxysize
 476               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
 477          enddo
 478 C Now Z dimension
 479 C Suggested by Ana -begins
 480         allareout=1
 481 C Suggested by Ana -ends
 482        j=3
 483         chain_beg=1
 484         do i=1,nres-1
 485          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 486           chain_end=i
 487           if (allareout.eq.1) then
 488             ireturnval=int(c(j,i)/boxysize)
 489             if (c(j,i).le.0) ireturnval=ireturnval-1
 490             do k=chain_beg,chain_end
 491               c(j,k)=c(j,k)-ireturnval*boxzsize
 492               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
 493             enddo
 494 C Suggested by Ana
 495             if (chain_beg.eq.1)
 496      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 497 C Suggested by Ana -end
 498            endif
 499            chain_beg=i+1
 500            allareout=1
 501          else
 502           if (int(c(j,i)/boxzsize).eq.0) allareout=0
 503          endif
 504         enddo
 505          if (allareout.eq.1) then
 506             ireturnval=int(c(j,i)/boxzsize)
 507             if (c(j,i).le.0) ireturnval=ireturnval-1
 508             do k=chain_beg,nres
 509               c(j,k)=c(j,k)-ireturnval*boxzsize
 510               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
 511             enddo
 512           endif
 513         nojumpval=0
 514         do i=2,nres
 515            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 516              difference=abs(c(j,i-1)-c(j,i))
 517              if (difference.gt.(boxzsize/2.0)) then
 518                 if (c(j,i-1).gt.c(j,i)) then
 519                   nojumpval=1
 520                  else
 521                    nojumpval=-1
 522                  endif
 523               else
 524               nojumpval=0
 525               endif
 526             endif
 527               c(j,i)=c(j,i)+nojumpval*boxzsize
 528               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
 529          enddo
 530        nojumpval=0
 531         do i=2,nres
 532            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 533              difference=abs(c(j,i-1)-c(j,i))
 534              if (difference.gt.boxzsize/2.0) then
 535                 if (c(j,i-1).gt.c(j,i)) then
 536                   nojumpval=1
 537                  else
 538                    nojumpval=-1
 539                  endif
 540               else
 541               nojumpval=0
 542               endif
 543             endif
 544              c(j,i)=c(j,i)+nojumpval*boxzsize
 545               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
 546          enddo
 547
 548         return
 549         end
 550