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 sc_coord_rebuild(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,rp(3),theta2,cost2,sint2,rj
 260       double precision scalar
 261       double precision ref(3,3),scalp,sscalp,refnorm
 262       double precision alpha,beta
 263 c      dsci=dsc(itype(i))
 264 c      dsci_inv=dsc_inv(itype(i))
 265       dsci=vbld(i+nres)
 266       dsci_inv=vbld_inv(i+nres)
 267 #ifdef OSF
 268       alphi=alph(i)
 269       omegi=omeg(i)
 270       if (alphi.ne.alphi) alphi=100.0
 271       if (omegi.ne.omegi) omegi=-100.0
 272 #else
 273       alphi=alph(i)
 274       omegi=omeg(i)
 275 #endif
 276       cosalphi=dcos(alphi)
 277       sinalphi=dsin(alphi)
 278       cosomegi=dcos(omegi)
 279       sinomegi=dsin(omegi)
 280       rp(1)= cosalphi
 281       rp(2)= sinalphi*cosomegi
 282       rp(3)=-sinalphi*sinomegi
 283 c Build the reference system
 284       do j=1,3
 285         ref(j,1)=-dc_norm(j,i-1)+dc_norm(j,i)
 286       enddo
 287       refnorm=dsqrt(scalar(ref(1,1),ref(1,1)))
 288       do j=1,3
 289         ref(j,1)=ref(j,1)/refnorm
 290       enddo
 291       scalp=scalar(ref(1,1),dc_norm(1,i))
 292       sscalp=1.0d0/dsqrt(1.0d0-scalp*scalp)
 293       do j=1,3
 294         ref(j,2)=(dc_norm(j,i)-scalp*ref(j,1))*sscalp
 295       enddo
 296       ref(1,3)= ref(2,1)*ref(3,2)-ref(3,1)*ref(2,2)
 297       ref(2,3)=-ref(1,1)*ref(3,2)+ref(3,1)*ref(1,2)
 298       ref(3,3)= ref(1,1)*ref(2,2)-ref(2,1)*ref(1,2)
 299 c      do j=1,3
 300 c        write (iout,*) j,scalar(ref(1,j),ref(1,1)),
 301 c     &  scalar(ref(1,j),ref(1,2)),scalar(ref(1,j),ref(1,3))
 302 c      enddo
 303 c Bring the coordinates to the global reference system
 304       do j=1,3
 305         dc_norm(j,nres+i)=0.0d0
 306         do k=1,3
 307           dc_norm(j,nres+i)=dc_norm(j,nres+i)+ref(j,k)*rp(k)
 308         enddo
 309         dc(j,nres+i)=dc_norm(j,nres+i)*dsci
 310         c(j,nres+i)=c(j,i)+dc(j,nres+i)
 311       enddo
 312 c      write (iout,*) scalar(dc_norm(1,i+nres),dc_norm(1,i+nres)),
 313 c     &  dsqrt(scalar(dc(1,i+nres),dc(1,i+nres)))
 314 c Check the internal coordinates
 315 c      c(:,2*nres+1)=ref(:,1)+c(:,i)
 316 c      write (iout,*) "alpha",rad2deg*alphi,
 317 c     & rad2deg*alpha(nres+i,i,2*nres+1)
 318 c      write (iout,*) "omega",rad2deg*omegi,
 319 c     & rad2deg*beta(nres+i,i,2*nres+1,i+1)
 320       return
 321       end
 322 c-----------------------------------------------------------------------------
 323       subroutine locate_side_chain(i)
 324 C
 325 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
 326 C
 327       implicit none
 328       include 'DIMENSIONS'
 329       include 'COMMON.CHAIN'
 330       include 'COMMON.LOCAL'
 331       include 'COMMON.GEO'
 332       include 'COMMON.VAR'
 333       include 'COMMON.IOUNITS'
 334       include 'COMMON.NAMES'
 335       include 'COMMON.INTERACT'
 336       integer i,j,k
 337       double precision xx(3)
 338       double precision dsci,dsci_inv,alphi,omegi,cosalphi,sinalphi,
 339      &  cosomegi,sinomegi,xp,yp,zp,theta2,cost2,sint2,rj
 340
 341 c      dsci=dsc(itype(i))
 342 c      dsci_inv=dsc_inv(itype(i))
 343       dsci=vbld(i+nres)
 344       dsci_inv=vbld_inv(i+nres)
 345 #ifdef OSF
 346       alphi=alph(i)
 347       omegi=omeg(i)
 348       if (alphi.ne.alphi) alphi=100.0
 349       if (omegi.ne.omegi) omegi=-100.0
 350 #else
 351       alphi=alph(i)
 352       omegi=omeg(i)
 353 #endif
 354       cosalphi=dcos(alphi)
 355       sinalphi=dsin(alphi)
 356       cosomegi=dcos(omegi)
 357       sinomegi=dsin(omegi)
 358       xp= dsci*cosalphi
 359       yp= dsci*sinalphi*cosomegi
 360       zp=-dsci*sinalphi*sinomegi
 361 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
 362 * X-axis aligned with the vector DC(*,i)
 363       theta2=pi-0.5D0*theta(i+1)
 364       cost2=dcos(theta2)
 365       sint2=dsin(theta2)
 366       xx(1)= xp*cost2+yp*sint2
 367       xx(2)=-xp*sint2+yp*cost2
 368       xx(3)= zp
 369 cd    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
 370 cd   &   xp,yp,zp,(xx(k),k=1,3)
 371       do j=1,3
 372         xloc(j,i)=xx(j)
 373       enddo
 374 * Bring the SC vectors to the common coordinate system.
 375       xx(1)=xloc(1,i)
 376       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
 377       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
 378       do j=1,3
 379         xrot(j,i)=xx(j)
 380       enddo
 381       do j=1,3
 382         rj=0.0D0
 383         do k=1,3
 384           rj=rj+prod(j,k,i-1)*xx(k)
 385         enddo
 386         dc(j,nres+i)=rj
 387         dc_norm(j,nres+i)=rj*dsci_inv
 388         c(j,nres+i)=c(j,i)+rj
 389       enddo
 390       return
 391       end
 392 c------------------------------------------
 393       subroutine returnbox
 394       include 'DIMENSIONS'
 395 #ifdef MPI
 396       include 'mpif.h'
 397 #endif
 398       include 'COMMON.CONTROL'
 399       include 'COMMON.VAR'
 400       include 'COMMON.MD'
 401 #ifndef LANG0
 402       include 'COMMON.LANGEVIN'
 403 #else
 404 #ifdef FIVEDIAG
 405       include 'COMMON.LANGEVIN.lang0.5diag'
 406 #else
 407       include 'COMMON.LANGEVIN.lang0'
 408 #endif
 409 #endif
 410       include 'COMMON.CHAIN'
 411       include 'COMMON.DERIV'
 412       include 'COMMON.GEO'
 413       include 'COMMON.LOCAL'
 414       include 'COMMON.INTERACT'
 415       include 'COMMON.IOUNITS'
 416       include 'COMMON.NAMES'
 417       include 'COMMON.TIME1'
 418       include 'COMMON.REMD'
 419       include 'COMMON.SETUP'
 420       include 'COMMON.MUCA'
 421       include 'COMMON.HAIRPIN'
 422 C change suggested by Ana - begin
 423       integer allareout
 424       integer i,j
 425 C change suggested by Ana - end
 426         j=1
 427         chain_beg=1
 428 C        do i=1,nres
 429 C       write(*,*) 'initial', i,j,c(j,i)
 430 C        enddo
 431 C change suggested by Ana - begin
 432         allareout=1
 433 C change suggested by Ana -end
 434         do i=1,nres-1
 435          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 436           chain_end=i
 437           if (allareout.eq.1) then
 438             ireturnval=int(c(j,i)/boxxsize)
 439             if (c(j,i).le.0) ireturnval=ireturnval-1
 440             do k=chain_beg,chain_end
 441               c(j,k)=c(j,k)-ireturnval*boxxsize
 442               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
 443             enddo
 444 C Suggested by Ana
 445             if (chain_beg.eq.1)
 446      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 447 C Suggested by Ana -end
 448            endif
 449            chain_beg=i+1
 450            allareout=1
 451          else
 452           if (int(c(j,i)/boxxsize).eq.0) allareout=0
 453          endif
 454         enddo
 455          if (allareout.eq.1) then
 456             ireturnval=int(c(j,i)/boxxsize)
 457             if (c(j,i).le.0) ireturnval=ireturnval-1
 458             do k=chain_beg,nres
 459               c(j,k)=c(j,k)-ireturnval*boxxsize
 460               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
 461             enddo
 462           endif
 463 C NO JUMP
 464 C        do i=1,nres
 465 C        write(*,*) 'befor no jump', i,j,c(j,i)
 466 C        enddo
 467         nojumpval=0
 468         do i=2,nres
 469            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 470              difference=abs(c(j,i-1)-c(j,i))
 471 C             print *,'diff', difference
 472              if (difference.gt.boxxsize/2.0) then
 473                 if (c(j,i-1).gt.c(j,i)) then
 474                   nojumpval=1
 475                  else
 476                    nojumpval=-1
 477                  endif
 478               else
 479               nojumpval=0
 480               endif
 481               endif
 482               c(j,i)=c(j,i)+nojumpval*boxxsize
 483               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
 484          enddo
 485        nojumpval=0
 486         do i=2,nres
 487            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 488              difference=abs(c(j,i-1)-c(j,i))
 489              if (difference.gt.boxxsize/2.0) then
 490                 if (c(j,i-1).gt.c(j,i)) then
 491                   nojumpval=1
 492                  else
 493                    nojumpval=-1
 494                  endif
 495               else
 496               nojumpval=0
 497               endif
 498              endif
 499               c(j,i)=c(j,i)+nojumpval*boxxsize
 500               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
 501          enddo
 502
 503 C        do i=1,nres
 504 C        write(*,*) 'after no jump', i,j,c(j,i)
 505 C        enddo
 506
 507 C NOW Y dimension
 508 C suggesed by Ana begins
 509         allareout=1
 510 C suggested by Ana ends
 511         j=2
 512         chain_beg=1
 513         do i=1,nres-1
 514          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 515           chain_end=i
 516           if (allareout.eq.1) then
 517             ireturnval=int(c(j,i)/boxysize)
 518             if (c(j,i).le.0) ireturnval=ireturnval-1
 519             do k=chain_beg,chain_end
 520               c(j,k)=c(j,k)-ireturnval*boxysize
 521              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
 522             enddo
 523 C Suggested by Ana
 524             if (chain_beg.eq.1)
 525      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 526 C Suggested by Ana -end
 527            endif
 528            chain_beg=i+1
 529            allareout=1
 530          else
 531           if (int(c(j,i)/boxysize).eq.0) allareout=0
 532          endif
 533         enddo
 534          if (allareout.eq.1) then
 535             ireturnval=int(c(j,i)/boxysize)
 536             if (c(j,i).le.0) ireturnval=ireturnval-1
 537             do k=chain_beg,nres
 538               c(j,k)=c(j,k)-ireturnval*boxysize
 539               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
 540             enddo
 541           endif
 542         nojumpval=0
 543         do i=2,nres
 544            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 545              difference=abs(c(j,i-1)-c(j,i))
 546              if (difference.gt.boxysize/2.0) then
 547                 if (c(j,i-1).gt.c(j,i)) then
 548                   nojumpval=1
 549                  else
 550                    nojumpval=-1
 551                  endif
 552              else
 553               nojumpval=0
 554               endif
 555            endif
 556               c(j,i)=c(j,i)+nojumpval*boxysize
 557               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
 558          enddo
 559       nojumpval=0
 560         do i=2,nres
 561            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 562              difference=abs(c(j,i-1)-c(j,i))
 563              if (difference.gt.boxysize/2.0) then
 564                 if (c(j,i-1).gt.c(j,i)) then
 565                   nojumpval=1
 566                  else
 567                    nojumpval=-1
 568                  endif
 569               else
 570               nojumpval=0
 571               endif
 572             endif
 573               c(j,i)=c(j,i)+nojumpval*boxysize
 574               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
 575          enddo
 576 C Now Z dimension
 577 C Suggested by Ana -begins
 578         allareout=1
 579 C Suggested by Ana -ends
 580        j=3
 581         chain_beg=1
 582         do i=1,nres-1
 583          if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
 584           chain_end=i
 585           if (allareout.eq.1) then
 586             ireturnval=int(c(j,i)/boxysize)
 587             if (c(j,i).le.0) ireturnval=ireturnval-1
 588             do k=chain_beg,chain_end
 589               c(j,k)=c(j,k)-ireturnval*boxzsize
 590               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
 591             enddo
 592 C Suggested by Ana
 593             if (chain_beg.eq.1)
 594      &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
 595 C Suggested by Ana -end
 596            endif
 597            chain_beg=i+1
 598            allareout=1
 599          else
 600           if (int(c(j,i)/boxzsize).eq.0) allareout=0
 601          endif
 602         enddo
 603          if (allareout.eq.1) then
 604             ireturnval=int(c(j,i)/boxzsize)
 605             if (c(j,i).le.0) ireturnval=ireturnval-1
 606             do k=chain_beg,nres
 607               c(j,k)=c(j,k)-ireturnval*boxzsize
 608               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
 609             enddo
 610           endif
 611         nojumpval=0
 612         do i=2,nres
 613            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 614              difference=abs(c(j,i-1)-c(j,i))
 615              if (difference.gt.(boxzsize/2.0)) then
 616                 if (c(j,i-1).gt.c(j,i)) then
 617                   nojumpval=1
 618                  else
 619                    nojumpval=-1
 620                  endif
 621               else
 622               nojumpval=0
 623               endif
 624             endif
 625               c(j,i)=c(j,i)+nojumpval*boxzsize
 626               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
 627          enddo
 628        nojumpval=0
 629         do i=2,nres
 630            if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
 631              difference=abs(c(j,i-1)-c(j,i))
 632              if (difference.gt.boxzsize/2.0) then
 633                 if (c(j,i-1).gt.c(j,i)) then
 634                   nojumpval=1
 635                  else
 636                    nojumpval=-1
 637                  endif
 638               else
 639               nojumpval=0
 640               endif
 641             endif
 642              c(j,i)=c(j,i)+nojumpval*boxzsize
 643               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
 644          enddo
 645
 646         return
 647         end
 648