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