source/unres/src-5hdiag-tmp/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 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 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 c      write (iout,*) "Calling chainbuild_extconf"
  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       implicit none
 303       include 'DIMENSIONS'
 304 #ifdef MPI
 305       include 'mpif.h'
 306 #endif
 307       include 'COMMON.CONTROL'
 308       include 'COMMON.VAR'
 309       include 'COMMON.MD'
 310 #ifndef LANG0
 311       include 'COMMON.LANGEVIN'
 312 #else
 313       include 'COMMON.LANGEVIN.lang0'
 314 #endif
 315       include 'COMMON.CHAIN'
 316       include 'COMMON.DERIV'
 317       include 'COMMON.GEO'
 318       include 'COMMON.LOCAL'
 319       include 'COMMON.INTERACT'
 320       include 'COMMON.IOUNITS'
 321       include 'COMMON.NAMES'
 322       include 'COMMON.TIME1'
 323       include 'COMMON.REMD'
 324       include 'COMMON.SETUP'
 325       include 'COMMON.MUCA'
 326       include 'COMMON.HAIRPIN'
 327 C change suggested by Ana - begin
 328       integer allareout,chain_beg,chain_end,i,j,k
 329       integer ireturnval,nojumpval
 330       double precision difference
 331 C change suggested by Ana - end
 332         j=1
 333         chain_beg=1
 334 C        do i=1,nres
 335 C       write(*,*) 'initial', i,j,c(j,i)
 336 C        enddo
 337 C change suggested by Ana - begin
 338         allareout=1
 339 C change suggested by Ana -end
 340         do i=1,nres-1
 341          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 342           chain_end=i
 343           if (allareout.eq.1) then
 344             ireturnval=int(c(j,i)/boxxsize)
 345             if (c(j,i).le.0) ireturnval=ireturnval-1
 346             do k=chain_beg,chain_end
 347               c(j,k)=c(j,k)-ireturnval*boxxsize
 348               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
 349             enddo
 350 C Suggested by Ana
 351             if (chain_beg.eq.1)
 352      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 353 C Suggested by Ana -end
 354            endif
 355            chain_beg=i+1
 356            allareout=1
 357          else
 358           if (int(c(j,i)/boxxsize).eq.0) allareout=0
 359          endif
 360         enddo
 361          if (allareout.eq.1) then
 362             ireturnval=int(c(j,i)/boxxsize)
 363             if (c(j,i).le.0) ireturnval=ireturnval-1
 364             do k=chain_beg,nres
 365               c(j,k)=c(j,k)-ireturnval*boxxsize
 366               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
 367             enddo
 368           endif
 369 C NO JUMP
 370 C        do i=1,nres
 371 C        write(*,*) 'befor no jump', i,j,c(j,i)
 372 C        enddo
 373         nojumpval=0
 374         do i=2,nres
 375            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 376              difference=abs(c(j,i-1)-c(j,i))
 377 C             print *,'diff', difference
 378              if (difference.gt.boxxsize/2.0) then
 379                 if (c(j,i-1).gt.c(j,i)) then
 380                   nojumpval=1
 381                  else
 382                    nojumpval=-1
 383                  endif
 384               else
 385               nojumpval=0
 386               endif
 387               endif
 388               c(j,i)=c(j,i)+nojumpval*boxxsize
 389               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
 390          enddo
 391        nojumpval=0
 392         do i=2,nres
 393            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 394              difference=abs(c(j,i-1)-c(j,i))
 395              if (difference.gt.boxxsize/2.0) then
 396                 if (c(j,i-1).gt.c(j,i)) then
 397                   nojumpval=1
 398                  else
 399                    nojumpval=-1
 400                  endif
 401               else
 402               nojumpval=0
 403               endif
 404              endif
 405               c(j,i)=c(j,i)+nojumpval*boxxsize
 406               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
 407          enddo
 408
 409 C        do i=1,nres
 410 C        write(*,*) 'after no jump', i,j,c(j,i)
 411 C        enddo
 412
 413 C NOW Y dimension
 414 C suggesed by Ana begins
 415         allareout=1
 416 C suggested by Ana ends
 417         j=2
 418         chain_beg=1
 419         do i=1,nres-1
 420          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 421           chain_end=i
 422           if (allareout.eq.1) then
 423             ireturnval=int(c(j,i)/boxysize)
 424             if (c(j,i).le.0) ireturnval=ireturnval-1
 425             do k=chain_beg,chain_end
 426               c(j,k)=c(j,k)-ireturnval*boxysize
 427              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
 428             enddo
 429 C Suggested by Ana
 430             if (chain_beg.eq.1)
 431      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 432 C Suggested by Ana -end
 433            endif
 434            chain_beg=i+1
 435            allareout=1
 436          else
 437           if (int(c(j,i)/boxysize).eq.0) allareout=0
 438          endif
 439         enddo
 440          if (allareout.eq.1) then
 441             ireturnval=int(c(j,i)/boxysize)
 442             if (c(j,i).le.0) ireturnval=ireturnval-1
 443             do k=chain_beg,nres
 444               c(j,k)=c(j,k)-ireturnval*boxysize
 445               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
 446             enddo
 447           endif
 448         nojumpval=0
 449         do i=2,nres
 450            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 451              difference=abs(c(j,i-1)-c(j,i))
 452              if (difference.gt.boxysize/2.0) then
 453                 if (c(j,i-1).gt.c(j,i)) then
 454                   nojumpval=1
 455                  else
 456                    nojumpval=-1
 457                  endif
 458              else
 459               nojumpval=0
 460               endif
 461            endif
 462               c(j,i)=c(j,i)+nojumpval*boxysize
 463               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
 464          enddo
 465       nojumpval=0
 466         do i=2,nres
 467            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 468              difference=abs(c(j,i-1)-c(j,i))
 469              if (difference.gt.boxysize/2.0) then
 470                 if (c(j,i-1).gt.c(j,i)) then
 471                   nojumpval=1
 472                  else
 473                    nojumpval=-1
 474                  endif
 475               else
 476               nojumpval=0
 477               endif
 478             endif
 479               c(j,i)=c(j,i)+nojumpval*boxysize
 480               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
 481          enddo
 482 C Now Z dimension
 483 C Suggested by Ana -begins
 484         allareout=1
 485 C Suggested by Ana -ends
 486        j=3
 487         chain_beg=1
 488         do i=1,nres-1
 489          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 490           chain_end=i
 491           if (allareout.eq.1) then
 492             ireturnval=int(c(j,i)/boxysize)
 493             if (c(j,i).le.0) ireturnval=ireturnval-1
 494             do k=chain_beg,chain_end
 495               c(j,k)=c(j,k)-ireturnval*boxzsize
 496               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
 497             enddo
 498 C Suggested by Ana
 499             if (chain_beg.eq.1)
 500      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 501 C Suggested by Ana -end
 502            endif
 503            chain_beg=i+1
 504            allareout=1
 505          else
 506           if (int(c(j,i)/boxzsize).eq.0) allareout=0
 507          endif
 508         enddo
 509          if (allareout.eq.1) then
 510             ireturnval=int(c(j,i)/boxzsize)
 511             if (c(j,i).le.0) ireturnval=ireturnval-1
 512             do k=chain_beg,nres
 513               c(j,k)=c(j,k)-ireturnval*boxzsize
 514               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
 515             enddo
 516           endif
 517         nojumpval=0
 518         do i=2,nres
 519            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 520              difference=abs(c(j,i-1)-c(j,i))
 521              if (difference.gt.(boxzsize/2.0)) then
 522                 if (c(j,i-1).gt.c(j,i)) then
 523                   nojumpval=1
 524                  else
 525                    nojumpval=-1
 526                  endif
 527               else
 528               nojumpval=0
 529               endif
 530             endif
 531               c(j,i)=c(j,i)+nojumpval*boxzsize
 532               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
 533          enddo
 534        nojumpval=0
 535         do i=2,nres
 536            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 537              difference=abs(c(j,i-1)-c(j,i))
 538              if (difference.gt.boxzsize/2.0) then
 539                 if (c(j,i-1).gt.c(j,i)) then
 540                   nojumpval=1
 541                  else
 542                    nojumpval=-1
 543                  endif
 544               else
 545               nojumpval=0
 546               endif
 547             endif
 548              c(j,i)=c(j,i)+nojumpval*boxzsize
 549               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
 550          enddo
 551
 552         return
 553         end
 554