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