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