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