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