source/unres/geometry.F90

   1                       module geometry
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use names
   5       use math
   6       use MPI_data
   7       use geometry_data
   8       use control_data
   9       use energy_data
  10       implicit none
  11 !-----------------------------------------------------------------------------
  12 ! commom.bounds
  13 !      common /bounds/
  14 !-----------------------------------------------------------------------------
  15 ! commom.chain
  16 !      common /chain/
  17 !      common /rotmat/
  18       real(kind=8),dimension(:,:,:),allocatable :: t,r !(3,3,maxres)
  19 !-----------------------------------------------------------------------------
  20 ! common.geo
  21 !      common /geo/
  22 !-----------------------------------------------------------------------------
  23 ! common.locmove
  24 !     Variables (set in init routine) never modified by local_move
  25 !      common /loc_const/
  26       integer :: init_called
  27       logical :: locmove_output
  28       real(kind=8) :: min_theta, max_theta
  29       real(kind=8) :: dmin2,dmax2
  30       real(kind=8) :: flag,small,small2
  31 !     Workspace for local_move
  32 !      common /loc_work/
  33       integer :: a_n,b_n,res_n
  34       real(kind=8),dimension(0:7) :: a_ang
  35       real(kind=8),dimension(0:3) :: b_ang
  36       real(kind=8),dimension(0:11) :: res_ang
  37       logical,dimension(0:2,0:7) :: a_tab
  38       logical,dimension(0:2,0:3) :: b_tab
  39       logical,dimension(0:2,0:2,0:11) :: res_tab
  40 !-----------------------------------------------------------------------------
  41 !      integer,dimension(:),allocatable :: itype_pdb !(maxres) initialize in molread
  42 !-----------------------------------------------------------------------------
  43 !
  44 !
  45 !-----------------------------------------------------------------------------
  46       contains
  47 !-----------------------------------------------------------------------------
  48 ! arcos.f
  49 !-----------------------------------------------------------------------------
  50       real(kind=8) function ARCOS(X)
  51 !      implicit real*8 (a-h,o-z)
  52 !      include 'COMMON.GEO'
  53 !el local variables
  54       real(kind=8) :: x
  55       IF (DABS(X).LT.1.0D0) GOTO 1
  56       ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
  57       RETURN
  58     1 ARCOS=DACOS(X)
  59       return
  60       end function ARCOS
  61 !-----------------------------------------------------------------------------
  62 ! chainbuild.F
  63 !-----------------------------------------------------------------------------
  64       subroutine chainbuild
  65 !
  66 ! Build the virtual polypeptide chain. Side-chain centroids are moveable.
  67 ! As of 2/17/95.
  68 !
  69 !      implicit real*8 (a-h,o-z)
  70 !      include 'DIMENSIONS'
  71 !      include 'COMMON.CHAIN'
  72 !      include 'COMMON.LOCAL'
  73 !      include 'COMMON.GEO'
  74 !      include 'COMMON.VAR'
  75 !      include 'COMMON.IOUNITS'
  76 !      include 'COMMON.NAMES'
  77 !      include 'COMMON.INTERACT'
  78       logical :: lprn
  79 !el local variables
  80       integer :: i,j
  81       real(kind=8) :: be,be1,alfai
  82       integer :: nres2
  83       nres2=2*nres
  84 ! Set lprn=.true. for debugging
  85       lprn = .false.
  86       print *,"I ENTER CHAINBUILD"
  87 !
  88 ! Define the origin and orientation of the coordinate system and locate the
  89 ! first three CA's and SC(2).
  90 !
  91 !elwrite(iout,*)"in chainbuild"
  92       call orig_frame
  93 !elwrite(iout,*)"after orig_frame"
  94 !
  95 ! Build the alpha-carbon chain.
  96 !
  97       do i=4,nres
  98         call locate_next_res(i)
  99       enddo
 100 !elwrite(iout,*)"after locate_next_res"
 101 !
 102 ! First and last SC must coincide with the corresponding CA.
 103 !
 104       do j=1,3
 105         dc(j,nres+1)=0.0D0
 106         dc_norm(j,nres+1)=0.0D0
 107         dc(j,nres+nres)=0.0D0
 108         dc_norm(j,nres+nres)=0.0D0
 109         c(j,nres+1)=c(j,1)
 110         c(j,nres+nres)=c(j,nres)
 111       enddo
 112 !
 113 ! Temporary diagnosis
 114 !
 115       if (lprn) then
 116
 117       call cartprint
 118       write (iout,'(/a)') 'Recalculated internal coordinates'
 119       do i=2,nres-1
 120         do j=1,3
 121           c(j,nres2+2)=0.5D0*(c(j,i-1)+c(j,i+1))        !maxres2=2*maxres
 122         enddo
 123         be=0.0D0
 124         if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
 125         be1=rad2deg*beta(nres+i,i,nres2+2,i+1)
 126         alfai=0.0D0
 127         if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
 128         write (iout,1212) restyp(itype(i,1),1),i,dist(i-1,i),&
 129         alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
 130       enddo
 131  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
 132
 133       endif
 134
 135       return
 136       end subroutine chainbuild
 137 !-----------------------------------------------------------------------------
 138       subroutine orig_frame
 139 !
 140 ! Define the origin and orientation of the coordinate system and locate
 141 ! the first three atoms.
 142 !
 143 !      implicit real*8 (a-h,o-z)
 144 !      include 'DIMENSIONS'
 145 !      include 'COMMON.CHAIN'
 146 !      include 'COMMON.LOCAL'
 147 !      include 'COMMON.GEO'
 148 !      include 'COMMON.VAR'
 149 !el local variables
 150       integer :: i,j
 151       real(kind=8) :: cost,sint
 152
 153 !el      allocate(t(3,3,nres))  !(3,3,maxres)
 154 !el      allocate(r(3,3,nres))  !(3,3,maxres)
 155 !el      allocate(rt(3,3,nres)) !(3,3,maxres)
 156 !el      allocate(dc_norm(3,0:2*nres))  !(3,0:maxres2)
 157 !el      allocate(prod(3,3,nres))       !(3,3,maxres)
 158
 159       cost=dcos(theta(3))
 160       sint=dsin(theta(3))
 161       t(1,1,1)=-cost
 162       t(1,2,1)=-sint
 163       t(1,3,1)= 0.0D0
 164       t(2,1,1)=-sint
 165       t(2,2,1)= cost
 166       t(2,3,1)= 0.0D0
 167       t(3,1,1)= 0.0D0
 168       t(3,2,1)= 0.0D0
 169       t(3,3,1)= 1.0D0
 170       r(1,1,1)= 1.0D0
 171       r(1,2,1)= 0.0D0
 172       r(1,3,1)= 0.0D0
 173       r(2,1,1)= 0.0D0
 174       r(2,2,1)= 1.0D0
 175       r(2,3,1)= 0.0D0
 176       r(3,1,1)= 0.0D0
 177       r(3,2,1)= 0.0D0
 178       r(3,3,1)= 1.0D0
 179       do i=1,3
 180         do j=1,3
 181           rt(i,j,1)=t(i,j,1)
 182         enddo
 183       enddo
 184       do i=1,3
 185         do j=1,3
 186           prod(i,j,1)=0.0D0
 187           prod(i,j,2)=t(i,j,1)
 188         enddo
 189         prod(i,i,1)=1.0D0
 190       enddo
 191       c(1,1)=0.0D0
 192       c(2,1)=0.0D0
 193       c(3,1)=0.0D0
 194       c(1,2)=vbld(2)
 195       c(2,2)=0.0D0
 196       c(3,2)=0.0D0
 197       dc(1,0)=0.0d0
 198       dc(2,0)=0.0D0
 199       dc(3,0)=0.0D0
 200       dc(1,1)=vbld(2)
 201       dc(2,1)=0.0D0
 202       dc(3,1)=0.0D0
 203       dc_norm(1,0)=0.0D0
 204       dc_norm(2,0)=0.0D0
 205       dc_norm(3,0)=0.0D0
 206       dc_norm(1,1)=1.0D0
 207       dc_norm(2,1)=0.0D0
 208       dc_norm(3,1)=0.0D0
 209       do j=1,3
 210         dc_norm(j,2)=prod(j,1,2)
 211         dc(j,2)=vbld(3)*prod(j,1,2)
 212         c(j,3)=c(j,2)+dc(j,2)
 213       enddo
 214       call locate_side_chain(2)
 215       return
 216       end subroutine orig_frame
 217 !-----------------------------------------------------------------------------
 218       subroutine locate_next_res(i)
 219 !
 220 ! Locate CA(i) and SC(i-1)
 221 !
 222 !      implicit real*8 (a-h,o-z)
 223 !      include 'DIMENSIONS'
 224 !      include 'COMMON.CHAIN'
 225 !      include 'COMMON.LOCAL'
 226 !      include 'COMMON.GEO'
 227 !      include 'COMMON.VAR'
 228 !      include 'COMMON.IOUNITS'
 229 !      include 'COMMON.NAMES'
 230 !      include 'COMMON.INTERACT'
 231 !
 232 ! Define the rotation matrices corresponding to CA(i)
 233 !
 234 !el local variables
 235       integer :: i,j
 236       real(kind=8) :: theti,phii
 237       real(kind=8) :: cost,sint,cosphi,sinphi
 238 #ifdef OSF
 239 #ifdef WHAM_RUN
 240       theti=theta(i)
 241       icrc=0
 242       call proc_proc(theti,icrc)
 243       if(icrc.eq.1)theti=100.0
 244       phii=phi(i)
 245       icrc=0
 246       call proc_proc(phii,icrc)
 247       if(icrc.eq.1)phii=180.0
 248 #else
 249       theti=theta(i)
 250       if (theti.ne.theti) theti=100.0
 251       phii=phi(i)
 252       if (phii.ne.phii) phii=180.0
 253 #endif
 254 #else
 255       theti=theta(i)
 256       phii=phi(i)
 257 #endif
 258       cost=dcos(theti)
 259       sint=dsin(theti)
 260       cosphi=dcos(phii)
 261       sinphi=dsin(phii)
 262 ! Define the matrices of the rotation about the virtual-bond valence angles
 263 ! theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
 264 ! program), R(i,j,k), and, the cumulative matrices of rotation RT
 265       t(1,1,i-2)=-cost
 266       t(1,2,i-2)=-sint
 267       t(1,3,i-2)= 0.0D0
 268       t(2,1,i-2)=-sint
 269       t(2,2,i-2)= cost
 270       t(2,3,i-2)= 0.0D0
 271       t(3,1,i-2)= 0.0D0
 272       t(3,2,i-2)= 0.0D0
 273       t(3,3,i-2)= 1.0D0
 274       r(1,1,i-2)= 1.0D0
 275       r(1,2,i-2)= 0.0D0
 276       r(1,3,i-2)= 0.0D0
 277       r(2,1,i-2)= 0.0D0
 278       r(2,2,i-2)=-cosphi
 279       r(2,3,i-2)= sinphi
 280       r(3,1,i-2)= 0.0D0
 281       r(3,2,i-2)= sinphi
 282       r(3,3,i-2)= cosphi
 283       rt(1,1,i-2)=-cost
 284       rt(1,2,i-2)=-sint
 285       rt(1,3,i-2)=0.0D0
 286       rt(2,1,i-2)=sint*cosphi
 287       rt(2,2,i-2)=-cost*cosphi
 288       rt(2,3,i-2)=sinphi
 289       rt(3,1,i-2)=-sint*sinphi
 290       rt(3,2,i-2)=cost*sinphi
 291       rt(3,3,i-2)=cosphi
 292       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
 293       do j=1,3
 294         dc_norm(j,i-1)=prod(j,1,i-1)
 295         dc(j,i-1)=vbld(i)*prod(j,1,i-1)
 296         c(j,i)=c(j,i-1)+dc(j,i-1)
 297       enddo
 298 !d    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
 299 !
 300 ! Now calculate the coordinates of SC(i-1)
 301 !
 302       call locate_side_chain(i-1)
 303       return
 304       end subroutine locate_next_res
 305 !-----------------------------------------------------------------------------
 306       subroutine locate_side_chain(i)
 307 !
 308 ! Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
 309 !
 310 !      implicit real*8 (a-h,o-z)
 311 !      include 'DIMENSIONS'
 312 !      include 'COMMON.CHAIN'
 313 !      include 'COMMON.LOCAL'
 314 !      include 'COMMON.GEO'
 315 !      include 'COMMON.VAR'
 316 !      include 'COMMON.IOUNITS'
 317 !      include 'COMMON.NAMES'
 318 !      include 'COMMON.INTERACT'
 319       integer :: i,j,k
 320       real(kind=8),dimension(3) :: xx
 321       real(kind=8) :: alphi,omegi,theta2
 322       real(kind=8) :: dsci,dsci_inv,sinalphi,cosalphi,cosomegi,sinomegi
 323       real(kind=8) :: xp,yp,zp,cost2,sint2,rj
 324 !      dsci=dsc(itype(i,1))
 325 !      dsci_inv=dsc_inv(itype(i,1))
 326       dsci=vbld(i+nres)
 327       dsci_inv=vbld_inv(i+nres)
 328 #ifdef OSF
 329       alphi=alph(i)
 330       omegi=omeg(i)
 331 #ifdef WHAM_RUN
 332 ! detecting NaNQ
 333       icrc=0
 334       call proc_proc(alphi,icrc)
 335       if(icrc.eq.1)alphi=100.0
 336       icrc=0
 337       call proc_proc(omegi,icrc)
 338       if(icrc.eq.1)omegi=-100.0
 339 #else
 340       if (alphi.ne.alphi) alphi=100.0
 341       if (omegi.ne.omegi) omegi=-100.0
 342 #endif
 343 #else
 344       alphi=alph(i)
 345       omegi=omeg(i)
 346 #endif
 347       cosalphi=dcos(alphi)
 348       sinalphi=dsin(alphi)
 349       cosomegi=dcos(omegi)
 350       sinomegi=dsin(omegi)
 351       xp= dsci*cosalphi
 352       yp= dsci*sinalphi*cosomegi
 353       zp=-dsci*sinalphi*sinomegi
 354 ! Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
 355 ! X-axis aligned with the vector DC(*,i)
 356       theta2=pi-0.5D0*theta(i+1)
 357       cost2=dcos(theta2)
 358       sint2=dsin(theta2)
 359       xx(1)= xp*cost2+yp*sint2
 360       xx(2)=-xp*sint2+yp*cost2
 361       xx(3)= zp
 362 !d    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i,1)),i,
 363 !d   &   xp,yp,zp,(xx(k),k=1,3)
 364       do j=1,3
 365         xloc(j,i)=xx(j)
 366       enddo
 367 ! Bring the SC vectors to the common coordinate system.
 368       xx(1)=xloc(1,i)
 369       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
 370       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
 371       do j=1,3
 372         xrot(j,i)=xx(j)
 373       enddo
 374       do j=1,3
 375         rj=0.0D0
 376         do k=1,3
 377           rj=rj+prod(j,k,i-1)*xx(k)
 378         enddo
 379         dc(j,nres+i)=rj
 380         dc_norm(j,nres+i)=rj*dsci_inv
 381         c(j,nres+i)=c(j,i)+rj
 382       enddo
 383       return
 384       end subroutine locate_side_chain
 385 !-----------------------------------------------------------------------------
 386 ! checkder_p.F
 387 !-----------------------------------------------------------------------------
 388       subroutine int_from_cart1(lprn)
 389 !      implicit real*8 (a-h,o-z)
 390 !      include 'DIMENSIONS'
 391 #ifdef MPI
 392       include 'mpif.h'
 393       integer :: ierror
 394 #endif
 395 !      include 'COMMON.IOUNITS'
 396 !      include 'COMMON.VAR'
 397 !      include 'COMMON.CHAIN'
 398 !      include 'COMMON.GEO'
 399 !      include 'COMMON.INTERACT'
 400 !      include 'COMMON.LOCAL'
 401 !      include 'COMMON.NAMES'
 402 !      include 'COMMON.SETUP'
 403 !      include 'COMMON.TIME1'
 404       logical :: lprn
 405 !el local variables
 406       integer :: i,j
 407       real(kind=8) :: dnorm1,dnorm2,be
 408       integer :: nres2
 409       nres2=2*nres
 410       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
 411 #ifdef TIMING
 412       time01=MPI_Wtime()
 413 #endif
 414
 415 #ifdef WHAM_RUN
 416       vbld(nres+1)=0.0d0
 417 !write(iout,*)"geometry warring, vbld=",(vbld(i),i=1,nres+1)
 418       vbld(2*nres)=0.0d0
 419       vbld_inv(nres+1)=0.0d0
 420       vbld_inv(2*nres)=0.0d0
 421 #endif
 422
 423 #if defined(PARINT) && defined(MPI)
 424       do i=iint_start,iint_end
 425 #else
 426       do i=2,nres
 427 #endif
 428         dnorm1=dist(i-1,i)
 429         dnorm2=dist(i,i+1)
 430         do j=1,3
 431           c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 &
 432            +(c(j,i+1)-c(j,i))/dnorm2)
 433         enddo
 434         be=0.0D0
 435         if (i.gt.2) then
 436         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
 437         if ((itype(i,1).ne.10).and.(itype(i-1,1).ne.10)) then
 438          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
 439         endif
 440         if (itype(i-1,1).ne.10) then
 441          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
 442          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
 443          omicron(2,i)=alpha(i-1+nres,i-1,i)
 444         endif
 445         if (itype(i,1).ne.10) then
 446          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
 447         endif
 448         endif
 449         omeg(i)=beta(nres+i,i,nres2+2,i+1)
 450         alph(i)=alpha(nres+i,i,nres2+2)
 451         theta(i+1)=alpha(i-1,i,i+1)
 452         vbld(i)=dist(i-1,i)
 453 !        print *,i,vbld(i),"vbld(i)"
 454         vbld_inv(i)=1.0d0/vbld(i)
 455         vbld(nres+i)=dist(nres+i,i)
 456         if ((itype(i,1).ne.10).and.(molnum(i).ne.5)) then
 457           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
 458         else
 459           vbld_inv(nres+i)=0.0d0
 460         endif
 461       enddo
 462 #if defined(PARINT) && defined(MPI)
 463        if (nfgtasks1.gt.1) then
 464 !d       write(iout,*) "iint_start",iint_start," iint_count",
 465 !d     &   (iint_count(i),i=0,nfgtasks-1)," iint_displ",
 466 !d     &   (iint_displ(i),i=0,nfgtasks-1)
 467 !d       write (iout,*) "Gather vbld backbone"
 468 !d       call flush(iout)
 469        time00=MPI_Wtime()
 470        call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),&
 471          MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),&
 472          MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 473 !d       write (iout,*) "Gather vbld_inv"
 474 !d       call flush(iout)
 475        call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),&
 476          MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),&
 477          MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 478 !d       write (iout,*) "Gather vbld side chain"
 479 !d       call flush(iout)
 480        call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),&
 481          MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),&
 482          MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 483 !d       write (iout,*) "Gather vbld_inv side chain"
 484 !d       call flush(iout)
 485        call MPI_Allgatherv(vbld_inv(iint_start+nres),&
 486          iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),&
 487          iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 488 !d       write (iout,*) "Gather theta"
 489 !d       call flush(iout)
 490        call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),&
 491          MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),&
 492          MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 493 !d       write (iout,*) "Gather phi"
 494 !d       call flush(iout)
 495        call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),&
 496          MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),&
 497          MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 498 #ifdef CRYST_SC
 499 !d       write (iout,*) "Gather alph"
 500 !d       call flush(iout)
 501        call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),&
 502          MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),&
 503          MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 504 !d       write (iout,*) "Gather omeg"
 505 !d       call flush(iout)
 506        call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),&
 507          MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),&
 508          MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 509 #endif
 510        time_gather=time_gather+MPI_Wtime()-time00
 511       endif
 512 #endif
 513       do i=1,nres-1
 514         do j=1,3
 515 !#ifdef WHAM_RUN
 516 #if defined(WHAM_RUN) || defined(CLUSTER)
 517           dc(j,i)=c(j,i+1)-c(j,i)
 518 #endif
 519           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 520         enddo
 521       enddo
 522       do i=2,nres-1
 523         do j=1,3
 524 !#ifdef WHAM_RUN
 525 #if defined(WHAM_RUN) || defined(CLUSTER)
 526           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 527 #endif
 528           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 529         enddo
 530       enddo
 531       if (lprn) then
 532       do i=2,nres
 533        write (iout,1212) restyp(itype(i,1),1),i,vbld(i),&
 534        rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),&
 535        rad2deg*alph(i),rad2deg*omeg(i)
 536       enddo
 537       endif
 538  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
 539 #ifdef TIMING
 540       time_intfcart=time_intfcart+MPI_Wtime()-time01
 541 #endif
 542       return
 543       end subroutine int_from_cart1
 544 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
 545 !-----------------------------------------------------------------------------
 546 ! check_sc_distr.f
 547 !-----------------------------------------------------------------------------
 548       subroutine check_sc_distr
 549 !      implicit real*8 (a-h,o-z)
 550 !      include 'DIMENSIONS'
 551 !      include 'COMMON.TIME1'
 552 !      include 'COMMON.INTERACT'
 553 !      include 'COMMON.NAMES'
 554 !      include 'COMMON.GEO'
 555 !      include 'COMMON.HEADER'
 556 !      include 'COMMON.CONTROL'
 557       logical :: fail
 558       real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
 559       real(kind=8) :: hrtime,mintime,sectime
 560       integer,parameter :: MaxSample=10000000
 561       real(kind=8),parameter :: delt=1.0D0/MaxSample
 562       real(kind=8),dimension(0:72,0:90) :: prob
 563 !el local variables
 564       integer :: it,i,j,isample,indal,indom
 565       real(kind=8) :: al,om,dV
 566       dV=2.0D0*5.0D0*deg2rad*deg2rad
 567       print *,'dv=',dv
 568       do 10 it=1,1
 569         if ((it.eq.10).or.(it.eq.ntyp1)) goto 10
 570         open (20,file=restyp(it,1)//'_distr.sdc',status='unknown')
 571         call gen_side(it,90.0D0 * deg2rad,al,om,fail,1)
 572         close (20)
 573         goto 10
 574         open (20,file=restyp(it,1)//'_distr1.sdc',status='unknown')
 575         do i=0,90
 576           do j=0,72
 577             prob(j,i)=0.0D0
 578           enddo
 579         enddo
 580         do isample=1,MaxSample
 581           call gen_side(it,90.0D0 * deg2rad,al,om,fail,1)
 582           indal=rad2deg*al/2
 583           indom=(rad2deg*om+180.0D0)/5
 584           prob(indom,indal)=prob(indom,indal)+delt
 585         enddo
 586         do i=45,90
 587           do j=0,72
 588             write (20,'(2f10.3,1pd15.5)') 2*i+0.0D0,5*j-180.0D0,&
 589                     prob(j,i)/dV
 590           enddo
 591         enddo
 592    10   continue
 593       return
 594       end subroutine check_sc_distr
 595 #endif
 596 !-----------------------------------------------------------------------------
 597 ! convert.f
 598 !-----------------------------------------------------------------------------
 599       subroutine geom_to_var(n,x)
 600 !
 601 ! Transfer the geometry parameters to the variable array.
 602 ! The positions of variables are as follows:
 603 ! 1. Virtual-bond torsional angles: 1 thru nres-3
 604 ! 2. Virtual-bond valence angles: nres-2 thru 2*nres-5
 605 ! 3. The polar angles alpha of local SC orientation: 2*nres-4 thru
 606 !    2*nres-4+nside
 607 ! 4. The torsional angles omega of SC orientation: 2*nres-4+nside+1
 608 !    thru 2*nre-4+2*nside
 609 !
 610 !      implicit real*8 (a-h,o-z)
 611 !      include 'DIMENSIONS'
 612 !      include 'COMMON.VAR'
 613 !      include 'COMMON.GEO'
 614 !      include 'COMMON.CHAIN'
 615       integer :: n,i
 616       real(kind=8),dimension(n) :: x
 617 !d    print *,'nres',nres,' nphi',nphi,' ntheta',ntheta,' nvar',nvar
 618       do i=4,nres
 619         x(i-3)=phi(i)
 620 !d      print *,i,i-3,phi(i)
 621       enddo
 622       if (n.eq.nphi) return
 623       do i=3,nres
 624         x(i-2+nphi)=theta(i)
 625 !d      print *,i,i-2+nphi,theta(i)
 626       enddo
 627       if (n.eq.nphi+ntheta) return
 628       do i=2,nres-1
 629         if (ialph(i,1).gt.0) then
 630           x(ialph(i,1))=alph(i)
 631           x(ialph(i,1)+nside)=omeg(i)
 632 !d        print *,i,ialph(i,1),ialph(i,1)+nside,alph(i),omeg(i)
 633         endif
 634       enddo
 635       return
 636       end subroutine geom_to_var
 637 !-----------------------------------------------------------------------------
 638       subroutine var_to_geom(n,x)
 639 !
 640 ! Update geometry parameters according to the variable array.
 641 !
 642 !      implicit real*8 (a-h,o-z)
 643 !      include 'DIMENSIONS'
 644 !      include 'COMMON.VAR'
 645 !      include 'COMMON.CHAIN'
 646 !      include 'COMMON.GEO'
 647 !      include 'COMMON.IOUNITS'
 648       integer :: n,i,ii
 649       real(kind=8),dimension(n) :: x
 650       logical :: change !,reduce
 651 !el      alph=0.0d0
 652 !el      omeg=0.0d0
 653 !el      phi=0.0d0
 654 !el      theta=0.0d0
 655
 656       change=reduce(x)
 657       if (n.gt.nphi+ntheta) then
 658         do i=1,nside
 659           ii=ialph(i,2)
 660           alph(ii)=x(nphi+ntheta+i)
 661           omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
 662 !elwrite(iout,*) "alph",ii,alph
 663 !elwrite(iout,*) "omeg",ii,omeg
 664         enddo
 665       endif
 666       do i=4,nres
 667         phi(i)=x(i-3)
 668 !elwrite(iout,*) "phi",i,phi
 669       enddo
 670       if (n.eq.nphi) return
 671       do i=3,nres
 672         theta(i)=x(i-2+nphi)
 673 !elwrite(iout,*) "theta",i,theta
 674         if (theta(i).eq.pi) theta(i)=0.99d0*pi
 675         x(i-2+nphi)=theta(i)
 676       enddo
 677       return
 678       end subroutine var_to_geom
 679 !-----------------------------------------------------------------------------
 680       logical function convert_side(alphi,omegi)
 681 !      implicit none
 682       real(kind=8) :: alphi,omegi
 683 !el      real(kind=8) :: pinorm
 684 !      include 'COMMON.GEO'
 685       convert_side=.false.
 686 ! Apply periodicity restrictions.
 687       if (alphi.gt.pi) then
 688         alphi=dwapi-alphi
 689         omegi=pinorm(omegi+pi)
 690         convert_side=.true.
 691       endif
 692       return
 693       end function convert_side
 694 !-----------------------------------------------------------------------------
 695       logical function reduce(x)
 696 !
 697 ! Apply periodic restrictions to variables.
 698 !
 699 !      implicit real*8 (a-h,o-z)
 700 !      include 'DIMENSIONS'
 701 !      include 'COMMON.VAR'
 702 !      include 'COMMON.CHAIN'
 703 !      include 'COMMON.GEO'
 704       logical :: zm,zmiana      !,convert_side
 705       real(kind=8),dimension(nvar) :: x
 706       integer :: i,ii,iii
 707       zmiana=.false.
 708       do i=4,nres
 709         x(i-3)=pinorm(x(i-3))
 710       enddo
 711       if (nvar.gt.nphi+ntheta) then
 712         do i=1,nside
 713           ii=nphi+ntheta+i
 714           iii=ii+nside
 715           x(ii)=thetnorm(x(ii))
 716           x(iii)=pinorm(x(iii))
 717 ! Apply periodic restrictions.
 718           zm=convert_side(x(ii),x(iii))
 719           zmiana=zmiana.or.zm
 720         enddo
 721       endif
 722       if (nvar.eq.nphi) return
 723       do i=3,nres
 724         ii=i-2+nphi
 725         iii=i-3
 726         x(ii)=dmod(x(ii),dwapi)
 727 ! Apply periodic restrictions.
 728         if (x(ii).gt.pi) then
 729           zmiana=.true.
 730           x(ii)=dwapi-x(ii)
 731           if (iii.gt.0) x(iii)=pinorm(x(iii)+pi)
 732           if (i.lt.nres) x(iii+1)=pinorm(x(iii+1)+pi)
 733           ii=ialph(i-1,1)
 734           if (ii.gt.0) then
 735             x(ii)=dmod(pi-x(ii),dwapi)
 736             x(ii+nside)=pinorm(-x(ii+nside))
 737             zm=convert_side(x(ii),x(ii+nside))
 738           endif
 739         else if (x(ii).lt.-pi) then
 740           zmiana=.true.
 741           x(ii)=dwapi+x(ii)
 742           ii=ialph(i-1,1)
 743           if (ii.gt.0) then
 744             x(ii)=dmod(pi-x(ii),dwapi)
 745             x(ii+nside)=pinorm(-pi-x(ii+nside))
 746             zm=convert_side(x(ii),x(ii+nside))
 747           endif
 748         else if (x(ii).lt.0.0d0) then
 749           zmiana=.true.
 750           x(ii)=-x(ii)
 751           if (iii.gt.0) x(iii)=pinorm(x(iii)+pi)
 752           if (i.lt.nres) x(iii+1)=pinorm(x(iii+1)+pi)
 753           ii=ialph(i-1,1)
 754           if (ii.gt.0) then
 755             x(ii+nside)=pinorm(-x(ii+nside))
 756             zm=convert_side(x(ii),x(ii+nside))
 757           endif
 758         endif
 759       enddo
 760       reduce=zmiana
 761       return
 762       end function reduce
 763 !-----------------------------------------------------------------------------
 764       real(kind=8) function thetnorm(x)
 765 ! This function puts x within [0,2Pi].
 766       implicit none
 767       real(kind=8) :: x,xx
 768 !      include 'COMMON.GEO'
 769       xx=dmod(x,dwapi)
 770       if (xx.lt.0.0d0) xx=xx+dwapi
 771       if (xx.gt.0.9999d0*pi) xx=0.9999d0*pi
 772       thetnorm=xx
 773       return
 774       end function thetnorm
 775 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
 776 !-----------------------------------------------------------------------------
 777       subroutine var_to_geom_restr(n,xx)
 778 !
 779 ! Update geometry parameters according to the variable array.
 780 !
 781 !      implicit real*8 (a-h,o-z)
 782 !      include 'DIMENSIONS'
 783 !      include 'COMMON.VAR'
 784 !      include 'COMMON.CHAIN'
 785 !      include 'COMMON.GEO'
 786 !      include 'COMMON.IOUNITS'
 787       integer :: n,i,ii
 788       real(kind=8),dimension(6*nres) :: x,xx !(maxvar) (maxvar=6*maxres)
 789       logical :: change !,reduce
 790
 791       call xx2x(x,xx)
 792       change=reduce(x)
 793       do i=1,nside
 794           ii=ialph(i,2)
 795           alph(ii)=x(nphi+ntheta+i)
 796           omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
 797       enddo
 798       do i=4,nres
 799         phi(i)=x(i-3)
 800       enddo
 801       do i=3,nres
 802         theta(i)=x(i-2+nphi)
 803         if (theta(i).eq.pi) theta(i)=0.99d0*pi
 804         x(i-2+nphi)=theta(i)
 805       enddo
 806       return
 807       end subroutine var_to_geom_restr
 808 !-----------------------------------------------------------------------------
 809 ! gen_rand_conf.F
 810 !-----------------------------------------------------------------------------
 811       subroutine gen_rand_conf(nstart,*)
 812 ! Generate random conformation or chain cut and regrowth.
 813       use mcm_data
 814       use random, only: iran_num,ran_number
 815 !      implicit real*8 (a-h,o-z)
 816 !      include 'DIMENSIONS'
 817 !      include 'COMMON.CHAIN'
 818 !      include 'COMMON.LOCAL'
 819 !      include 'COMMON.VAR'
 820 !      include 'COMMON.INTERACT'
 821 !      include 'COMMON.IOUNITS'
 822 !      include 'COMMON.MCM'
 823 !      include 'COMMON.GEO'
 824 !      include 'COMMON.CONTROL'
 825       logical :: back,fail      !overlap,
 826 !el local variables
 827       integer :: i,nstart,maxsi,nsi,maxnit,nit,niter
 828       integer :: it1,it2,it,j
 829 !d    print *,' CG Processor',me,' maxgen=',maxgen
 830       maxsi=1000
 831       write (iout,*) 'Gen_Rand_conf: nstart=',nstart,nres
 832       if (nstart.lt.5) then
 833         it1=iabs(itype(2,1))
 834         phi(4)=gen_phi(4,iabs(itype(2,1)),iabs(itype(3,1)))
 835 !       write(iout,*)'phi(4)=',rad2deg*phi(4)
 836         if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4),molnum(2))
 837 !       write(iout,*)'theta(3)=',rad2deg*theta(3)
 838         if ((it1.ne.10).and.(it1.ne.ntyp1)) then
 839           nsi=0
 840           fail=.true.
 841           do while (fail.and.nsi.le.maxsi)
 842             call gen_side(it1,theta(3),alph(2),omeg(2),fail,molnum(2))
 843             write (iout,*) 'nsi=',nsi,maxsi
 844             nsi=nsi+1
 845           enddo
 846           if (nsi.gt.maxsi) return 1
 847         endif ! it1.ne.10
 848         write(iout,*) "before origin_frame"
 849         call orig_frame
 850         write(iout,*) "after origin_frame"
 851         i=4
 852         nstart=4
 853       else
 854         i=nstart
 855         nstart=max0(i,4)
 856       endif
 857
 858       maxnit=0
 859
 860       nit=0
 861       niter=0
 862       back=.false.
 863       do while (i.le.nres .and. niter.lt.maxgen)
 864         write(iout,*) 'i=',i,'back=',back
 865         if (i.lt.nstart) then
 866           if(iprint.gt.1) then
 867           write (iout,'(/80(1h*)/2a/80(1h*))') &
 868                 'Generation procedure went down to ',&
 869                 'chain beginning. Cannot continue...'
 870           write (*,'(/80(1h*)/2a/80(1h*))') &
 871                 'Generation procedure went down to ',&
 872                 'chain beginning. Cannot continue...'
 873           endif
 874           return 1
 875         endif
 876         it1=iabs(itype(i-1,molnum(i-1)))
 877         it2=iabs(itype(i-2,molnum(i-2)))
 878         it=iabs(itype(i,molnum(i)))
 879         if ((it.eq.ntyp1).and.(it1.eq.ntyp1)) &
 880           vbld(i)=ran_number(30.0D0,40.0D0)
 881 !       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,&
 882 !        ' nit=',nit,' niter=',niter,' maxgen=',maxgen
 883         phi(i+1)=gen_phi(i+1,it1,it)
 884         if (back) then
 885           phi(i)=gen_phi(i+1,it2,it1)
 886 !         print *,'phi(',i,')=',phi(i)
 887           theta(i-1)=gen_theta(it2,phi(i-1),phi(i),molnum(i))
 888 !          print *,"theta",theta(i-1),phi(i)
 889           if ((it2.ne.10).and.(it2.ne.ntyp1)) then
 890             nsi=0
 891             fail=.true.
 892             do while (fail.and.nsi.le.maxsi)
 893               call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail,molnum(i-2))
 894               nsi=nsi+1
 895             enddo
 896             if (nsi.gt.maxsi) return 1
 897           endif
 898           call locate_next_res(i-1)
 899         endif
 900         theta(i)=gen_theta(it1,phi(i),phi(i+1),molnum(i))
 901 !        write(iout,*) "theta(i),",theta(i)
 902         if ((it1.ne.10).and.(it1.ne.ntyp1)) then
 903         nsi=0
 904         fail=.true.
 905         do while (fail.and.nsi.le.maxsi)
 906           call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail,molnum(i))
 907 !                  write(iout,*)"alpha,omeg(i-1)",alph(i-1),omeg(i-1),i,nsi,maxsi
 908           nsi=nsi+1
 909         enddo
 910         if (nsi.gt.maxsi) return 1
 911         endif
 912         call locate_next_res(i)
 913         write(iout,*) "overlap,",overlap(i-1)
 914         if (overlap(i-1)) then
 915           if (nit.lt.maxnit) then
 916             back=.true.
 917             nit=nit+1
 918           else
 919             nit=0
 920             if (i.gt.3) then
 921               back=.true.
 922               i=i-1
 923             else
 924               write (iout,'(a)') &
 925         'Cannot generate non-overlaping conformation. Increase MAXNIT.'
 926               write (*,'(a)') &
 927         'Cannot generate non-overlaping conformation. Increase MAXNIT.'
 928               return 1
 929             endif
 930           endif
 931         else
 932 !          write(iout,*) "tu dochodze"
 933           back=.false.
 934           nit=0
 935           i=i+1
 936         endif
 937         niter=niter+1
 938       enddo
 939       if (niter.ge.maxgen) then
 940         write (iout,'(a,2i5)') &
 941        'Too many trials in conformation generation',niter,maxgen
 942         write (*,'(a,2i5)') &
 943        'Too many trials in conformation generation',niter,maxgen
 944         return 1
 945       endif
 946       do j=1,3
 947         c(j,nres+1)=c(j,1)
 948         c(j,nres+nres)=c(j,nres)
 949       enddo
 950       return
 951       end subroutine gen_rand_conf
 952 !-----------------------------------------------------------------------------
 953       logical function overlap(i)
 954 !      implicit real*8 (a-h,o-z)
 955 !      include 'DIMENSIONS'
 956 !      include 'COMMON.CHAIN'
 957 !      include 'COMMON.INTERACT'
 958 !      include 'COMMON.FFIELD'
 959       integer :: i,j,iti,itj,iteli,itelj,k
 960       real(kind=8) :: redfac,rcomp
 961       integer :: nres2
 962       nres2=2*nres
 963       data redfac /0.5D0/
 964       overlap=.false.
 965       iti=iabs(itype(i,molnum(i)))
 966       if (iti.gt.ntyp) return
 967 ! Check for SC-SC overlaps.
 968 !d    print *,'nnt=',nnt,' nct=',nct
 969       do j=nnt,i-1
 970 !        print *, "molnum(j)",j,molnum(j)
 971         if (molnum(j).eq.1) then
 972         itj=iabs(itype(j,1))
 973         if (itj.eq.ntyp1) cycle
 974         if (j.lt.i-1 .or. ipot.ne.4) then
 975           rcomp=sigmaii(iti,itj)
 976         else
 977           rcomp=sigma(iti,itj)
 978         endif
 979 !d      print *,'j=',j
 980         if (dist(nres+i,nres+j).lt.redfac*rcomp) then
 981           overlap=.true.
 982
 983 !        print *,'overlap, SC-SC: i=',i,' j=',j,
 984 !     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
 985 !     &     rcomp
 986           return
 987         endif
 988         else if (molnum(j).eq.2) then
 989         itj=iabs(itype(j,2))
 990         if (dist(nres+i,nres+j).lt.redfac*sigma_nucl(iti,itj)) then
 991           overlap=.true.
 992
 993 !        print *,'overlap, SC-SC: i=',i,' j=',j,
 994 !     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
 995 !     &     rcomp
 996           return
 997         endif
 998
 999       endif
1000       enddo
1001 ! Check for overlaps between the added peptide group and the preceding
1002 ! SCs.
1003       iteli=itel(i)
1004       do j=1,3
1005 !       c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1))
1006         c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
1007       enddo
1008       do j=nnt,i-2
1009         if (molnum(j).ne.1) cycle
1010         itj=iabs(itype(j,1))
1011 !d      print *,'overlap, p-Sc: i=',i,' j=',j,
1012 !d   &         ' dist=',dist(nres+j,maxres2+1)
1013         if (dist(nres+j,nres2+3).lt.4.0D0*redfac) then
1014           overlap=.true.
1015           return
1016         endif
1017       enddo
1018 ! Check for overlaps between the added side chain and the preceding peptide
1019 ! groups.
1020       do j=1,nnt-2
1021         if (molnum(j).ne.1) cycle
1022         do k=1,3
1023           c(k,nres2+3)=0.5D0*(c(k,j)+c(k,j+1))
1024         enddo
1025 !d      print *,'overlap, SC-p: i=',i,' j=',j,
1026 !d   &         ' dist=',dist(nres+i,maxres2+1)
1027         if (dist(nres+i,nres2+3).lt.4.0D0*redfac) then
1028           overlap=.true.
1029           return
1030         endif
1031       enddo
1032 ! Check for p-p overlaps
1033       do j=1,3
1034         c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
1035       enddo
1036       do j=nnt,i-2
1037 !        if (molnum(j).eq.1) then
1038         itelj=itel(j)
1039         do k=1,3
1040           c(k,nres2+4)=0.5D0*(c(k,j)+c(k,j+1))
1041         enddo
1042 !d      print *,'overlap, p-p: i=',i,' j=',j,
1043 !d   &         ' dist=',dist(maxres2+1,maxres2+2)
1044         if (molnum(j).eq.1) then
1045         if(iteli.ne.0.and.itelj.ne.0)then
1046         if (dist(nres2+3,nres2+4).lt.rpp(iteli,itelj)*redfac) then
1047           overlap=.true.
1048           return
1049         endif
1050         endif
1051         else if (molnum(j).eq.2) then
1052         if (dist(nres2+3,nres2+4).lt.3.0) then
1053           overlap=.true.
1054           return
1055         endif
1056       endif
1057       enddo
1058       return
1059       end function overlap
1060 !-----------------------------------------------------------------------------
1061       real(kind=8) function gen_phi(i,it1,it2)
1062       use random, only:ran_number
1063 !      implicit real*8 (a-h,o-z)
1064 !      include 'DIMENSIONS'
1065 !      include 'COMMON.GEO'
1066 !      include 'COMMON.BOUNDS'
1067       integer :: i,it1,it2
1068 !      gen_phi=ran_number(-pi,pi)
1069 ! 8/13/98 Generate phi using pre-defined boundaries
1070       gen_phi=ran_number(phibound(1,i),phibound(2,i))
1071       return
1072       end function gen_phi
1073 !-----------------------------------------------------------------------------
1074       real(kind=8) function gen_theta(it,gama,gama1,mnum)
1075       use random,only:binorm,ran_number
1076 !      implicit real*8 (a-h,o-z)
1077 !      include 'DIMENSIONS'
1078 !      include 'COMMON.LOCAL'
1079 !      include 'COMMON.GEO'
1080       real(kind=8),dimension(2) :: y,z
1081       real(kind=8) :: theta_max,theta_min,sig,ak
1082 !el local variables
1083       integer :: j,it,k,mnum
1084       real(kind=8) :: gama,gama1,thet_pred_mean,theta_temp
1085 !     print *,'gen_theta: it=',it
1086       theta_min=0.05D0*pi
1087       theta_max=0.95D0*pi
1088       if (dabs(gama).gt.dwapi) then
1089         y(1)=dcos(gama)
1090         y(2)=dsin(gama)
1091       else
1092         y(1)=0.0D0
1093         y(2)=0.0D0
1094       endif
1095       if (dabs(gama1).gt.dwapi) then
1096         z(1)=dcos(gama1)
1097         z(2)=dsin(gama1)
1098       else
1099         z(1)=0.0D0
1100         z(2)=0.0D0
1101       endif
1102       if (it.eq.ntyp1) then
1103       gen_theta=ran_number(theta_max/2.0,theta_max)
1104       else if (mnum.eq.1) then
1105
1106       thet_pred_mean=a0thet(it)
1107 !      write(iout,*),it,thet_pred_mean,"gen_thet"
1108       do k=1,2
1109         thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k) &
1110            +bthet(k,it,1,1)*z(k)
1111       enddo
1112       sig=polthet(3,it)
1113       do j=2,0,-1
1114         sig=sig*thet_pred_mean+polthet(j,it)
1115       enddo
1116       sig=0.5D0/(sig*sig+sigc0(it))
1117       ak=dexp(gthet(1,it)- &
1118        0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
1119 !     print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
1120 !     print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
1121       theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak)
1122       if (theta_temp.lt.theta_min) theta_temp=theta_min
1123       if (theta_temp.gt.theta_max) theta_temp=theta_max
1124       gen_theta=theta_temp
1125 !     print '(a)','Exiting GENTHETA.'
1126       else if (mnum.eq.2) then
1127        gen_theta=2.0d0 + ran_number(0.0d0,0.34d0)
1128       else
1129               gen_theta=ran_number(theta_max/2.0,theta_max)
1130        endif
1131       return
1132       end function gen_theta
1133 !-----------------------------------------------------------------------------
1134       subroutine gen_side(it,the,al,om,fail,mnum)
1135       use random, only:ran_number,mult_norm1
1136 !      implicit real*8 (a-h,o-z)
1137 !      include 'DIMENSIONS'
1138 !      include 'COMMON.GEO'
1139 !      include 'COMMON.LOCAL'
1140 !      include 'COMMON.SETUP'
1141 !      include 'COMMON.IOUNITS'
1142       real(kind=8) :: MaxBoxLen=10.0D0
1143       real(kind=8),dimension(3,3) :: Ap_inv,a,vec
1144       real(kind=8),dimension(:,:),allocatable :: z !(3,maxlob)
1145       real(kind=8),dimension(:),allocatable :: W1,detAp !(maxlob)
1146       real(kind=8),dimension(:),allocatable :: sumW !(0:maxlob)
1147       real(kind=8),dimension(2) :: y,cm,eig
1148       real(kind=8),dimension(2,2) :: box
1149       real(kind=8),dimension(100) :: work
1150       real(kind=8) :: eig_limit=1.0D-8
1151       real(kind=8) :: Big=10.0D0
1152       logical :: lprint,fail,lcheck
1153 !el local variables
1154       integer :: it,i,j,k,l,nlobit,ial,iom,iii,ilob,mnum
1155       real(kind=8) :: the,al,om,detApi,wart,y2,wykl,radmax
1156       real(kind=8) :: tant,zz1,W1i,radius,zk,fac,dV,sum,sum1
1157       real(kind=8) :: which_lobe
1158       lcheck=.false.
1159       lprint=.false.
1160       fail=.false.
1161       if (mnum.eq.1) then
1162       if (the.eq.0.0D0 .or. the.eq.pi) then
1163 #ifdef MPI
1164         write (*,'(a,i4,a,i3,a,1pe14.5)') &
1165        'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
1166 #else
1167 !d        write (iout,'(a,i3,a,1pe14.5)')
1168 !d     &   'Error in GenSide: it=',it,' theta=',the
1169 #endif
1170         fail=.true.
1171         return
1172       endif
1173       if (nlobit.eq.0) then
1174          al=ran_number(0.05d0,pi/2)
1175          om=ran_number(-pi,pi)
1176          return
1177       endif
1178       tant=dtan(the-pipol)
1179       nlobit=nlob(it)
1180       allocate(z(3,nlobit))
1181       allocate(W1(nlobit))
1182       allocate(detAp(nlobit))
1183       allocate(sumW(0:nlobit))
1184       if (lprint) then
1185 #ifdef MPI
1186         print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
1187         write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
1188 #endif
1189         print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
1190         write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,&
1191            ' tant=',tant
1192       endif
1193       do i=1,nlobit
1194        zz1=tant-censc(1,i,it)
1195         do k=1,3
1196           do l=1,3
1197             a(k,l)=gaussc(k,l,i,it)
1198           enddo
1199         enddo
1200         detApi=a(2,2)*a(3,3)-a(2,3)**2
1201         Ap_inv(2,2)=a(3,3)/detApi
1202         Ap_inv(2,3)=-a(2,3)/detApi
1203         Ap_inv(3,2)=Ap_inv(2,3)
1204         Ap_inv(3,3)=a(2,2)/detApi
1205         if (lprint) then
1206           write (*,'(/a,i2/)') 'Cluster #',i
1207           write (*,'(3(1pe14.5),5x,1pe14.5)') &
1208           ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
1209           write (iout,'(/a,i2/)') 'Cluster #',i
1210           write (iout,'(3(1pe14.5),5x,1pe14.5)') &
1211           ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
1212         endif
1213         W1i=0.0D0
1214         do k=2,3
1215           do l=2,3
1216             W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
1217           enddo
1218         enddo
1219         W1i=a(1,1)-W1i
1220         W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
1221 !        if (lprint) write(*,'(a,3(1pe15.5)/)')
1222 !     &          'detAp, W1, anormi',detApi,W1i,anormi
1223         do k=2,3
1224           zk=censc(k,i,it)
1225           do l=2,3
1226             zk=zk+zz1*Ap_inv(k,l)*a(l,1)
1227           enddo
1228           z(k,i)=zk
1229         enddo
1230         detAp(i)=dsqrt(detApi)
1231       enddo
1232
1233       if (lprint) then
1234         print *,'W1:',(w1(i),i=1,nlobit)
1235         print *,'detAp:',(detAp(i),i=1,nlobit)
1236         print *,'Z'
1237         do i=1,nlobit
1238           print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
1239         enddo
1240         write (iout,*) 'W1:',(w1(i),i=1,nlobit)
1241         write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
1242         write (iout,*) 'Z'
1243         do i=1,nlobit
1244           write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
1245         enddo
1246       endif
1247       if (lcheck) then
1248 ! Writing the distribution just to check the procedure
1249       fac=0.0D0
1250       dV=deg2rad**2*10.0D0
1251       sum=0.0D0
1252       sum1=0.0D0
1253       do i=1,nlobit
1254         fac=fac+W1(i)/detAp(i)
1255       enddo
1256       fac=1.0D0/(2.0D0*fac*pi)
1257 !d    print *,it,'fac=',fac
1258       do ial=90,180,2
1259         y(1)=deg2rad*ial
1260         do iom=-180,180,5
1261           y(2)=deg2rad*iom
1262           wart=0.0D0
1263           do i=1,nlobit
1264             do j=2,3
1265               do k=2,3
1266                 a(j-1,k-1)=gaussc(j,k,i,it)
1267               enddo
1268             enddo
1269             y2=y(2)
1270
1271             do iii=-1,1
1272
1273               y(2)=y2+iii*dwapi
1274
1275               wykl=0.0D0
1276               do j=1,2
1277                 do k=1,2
1278                   wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
1279                 enddo
1280               enddo
1281               wart=wart+W1(i)*dexp(-0.5D0*wykl)
1282
1283             enddo
1284
1285             y(2)=y2
1286
1287           enddo
1288 !         print *,'y',y(1),y(2),' fac=',fac
1289           wart=fac*wart
1290           write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
1291           sum=sum+wart
1292           sum1=sum1+1.0D0
1293         enddo
1294       enddo
1295 !     print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
1296       return
1297       endif
1298
1299 ! Calculate the CM of the system
1300 !
1301       do i=1,nlobit
1302         W1(i)=W1(i)/detAp(i)
1303       enddo
1304       sumW(0)=0.0D0
1305       do i=1,nlobit
1306         sumW(i)=sumW(i-1)+W1(i)
1307       enddo
1308       cm(1)=z(2,1)*W1(1)
1309       cm(2)=z(3,1)*W1(1)
1310       do j=2,nlobit
1311         cm(1)=cm(1)+z(2,j)*W1(j)
1312         cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
1313       enddo
1314       cm(1)=cm(1)/sumW(nlobit)
1315       cm(2)=cm(2)/sumW(nlobit)
1316       if (cm(1).gt.Big .or. cm(1).lt.-Big .or. &
1317        cm(2).gt.Big .or. cm(2).lt.-Big) then
1318 !d        write (iout,'(a)')
1319 !d     & 'Unexpected error in GenSide - CM coordinates too large.'
1320 !d        write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
1321 !d        write (*,'(a)')
1322 !d     & 'Unexpected error in GenSide - CM coordinates too large.'
1323 !d        write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
1324         fail=.true.
1325         return
1326       endif
1327 !d    print *,'CM:',cm(1),cm(2)
1328 !
1329 ! Find the largest search distance from CM
1330 !
1331       radmax=0.0D0
1332       do i=1,nlobit
1333         do j=2,3
1334           do k=2,3
1335             a(j-1,k-1)=gaussc(j,k,i,it)
1336           enddo
1337         enddo
1338 #ifdef NAG
1339         call f02faf('N','U',2,a,3,eig,work,100,ifail)
1340 #else
1341         call djacob(2,3,10000,1.0d-10,a,vec,eig)
1342 #endif
1343 #ifdef MPI
1344         if (lprint) then
1345           print *,'*************** CG Processor',me
1346           print *,'CM:',cm(1),cm(2)
1347           write (iout,*) '*************** CG Processor',me
1348           write (iout,*) 'CM:',cm(1),cm(2)
1349           print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
1350           write (iout,'(A,8f10.5)') &
1351               'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
1352         endif
1353 #endif
1354         if (eig(1).lt.eig_limit) then
1355           write(iout,'(a)') &
1356            'From Mult_Norm: Eigenvalues of A are too small.'
1357           write(*,'(a)') &
1358            'From Mult_Norm: Eigenvalues of A are too small.'
1359           fail=.true.
1360           return
1361         endif
1362         radius=0.0D0
1363 !d      print *,'i=',i
1364         do j=1,2
1365           radius=radius+pinorm(z(j+1,i)-cm(j))**2
1366         enddo
1367         radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
1368         if (radius.gt.radmax) radmax=radius
1369       enddo
1370       if (radmax.gt.pi) radmax=pi
1371 !
1372 ! Determine the boundaries of the search rectangle.
1373 !
1374       if (lprint) then
1375         print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
1376         print '(a,4(1pe14.4))','radmax: ',radmax
1377       endif
1378       box(1,1)=dmax1(cm(1)-radmax,0.0D0)
1379       box(2,1)=dmin1(cm(1)+radmax,pi)
1380       box(1,2)=cm(2)-radmax
1381       box(2,2)=cm(2)+radmax
1382       if (lprint) then
1383 #ifdef MPI
1384         print *,'CG Processor',me,' Array BOX:'
1385 #else
1386         print *,'Array BOX:'
1387 #endif
1388         print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
1389         print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
1390 #ifdef MPI
1391         write (iout,*)'CG Processor',me,' Array BOX:'
1392 #else
1393         write (iout,*)'Array BOX:'
1394 #endif
1395         write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
1396         write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
1397       endif
1398 !      if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
1399 !#ifdef MPI
1400 !        write (iout,'(a,i4,a,3e15.5)') 'CG Processor:',me,': bad sampling box.',box(1,2),box(2,2),radmax
1401 !        write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
1402 !#else
1403 !        write (iout,'(a)') 'Bad sampling box.'
1404 !#endif
1405 !        fail=.true.
1406 !        return
1407 !      endif
1408       which_lobe=ran_number(0.0D0,sumW(nlobit))
1409 !     print '(a,1pe14.4)','which_lobe=',which_lobe
1410       do i=1,nlobit
1411         if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
1412       enddo
1413     1 ilob=i
1414 !     print *,'ilob=',ilob,' nlob=',nlob(it)
1415       do i=2,3
1416         cm(i-1)=z(i,ilob)
1417         do j=2,3
1418           a(i-1,j-1)=gaussc(i,j,ilob,it)
1419         enddo
1420       enddo
1421 !d    print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
1422       call mult_norm1(3,2,a,cm,box,y,fail)
1423       if (fail) return
1424       al=y(1)
1425       om=pinorm(y(2))
1426       else if (mnum.eq.2) then
1427        al=0.7+ran_number(0.0d0,0.2d0)
1428        om=ran_number(0.0d0,3.14d0)
1429       endif
1430
1431 !d    print *,'al=',al,' om=',om
1432 !d    stop
1433       return
1434       end subroutine gen_side
1435 !-----------------------------------------------------------------------------
1436       subroutine overlap_sc(scfail)
1437 !
1438 !     Internal and cartesian coordinates must be consistent as input,
1439 !     and will be up-to-date on return.
1440 !     At the end of this procedure, scfail is true if there are
1441 !     overlapping residues left, or false otherwise (success)
1442 !
1443 !      implicit real*8 (a-h,o-z)
1444 !      include 'DIMENSIONS'
1445 !      include 'COMMON.CHAIN'
1446 !      include 'COMMON.INTERACT'
1447 !      include 'COMMON.FFIELD'
1448 !      include 'COMMON.VAR'
1449 !      include 'COMMON.SBRIDGE'
1450 !      include 'COMMON.IOUNITS'
1451       logical :: had_overlaps,fail,scfail
1452       integer,dimension(nres) :: ioverlap !(maxres)
1453       integer :: ioverlap_last,k,maxsi,i,iti,nsi
1454       integer :: ires,j
1455
1456       had_overlaps=.false.
1457       call overlap_sc_list(ioverlap,ioverlap_last)
1458       if (ioverlap_last.gt.0) then
1459         write (iout,*) '#OVERLAPing residues ',ioverlap_last
1460         write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
1461         had_overlaps=.true.
1462       endif
1463
1464       maxsi=1000
1465       do k=1,1000
1466         if (ioverlap_last.eq.0) exit
1467
1468         do ires=1,ioverlap_last
1469           i=ioverlap(ires)
1470           iti=iabs(itype(i,1))
1471           if ((iti.ne.10).and.(molnum(i).ne.5).and.(iti.ne.ntyp1)) then
1472             nsi=0
1473             fail=.true.
1474             do while (fail.and.nsi.le.maxsi)
1475               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i))
1476               nsi=nsi+1
1477             enddo
1478             if(fail) goto 999
1479           endif
1480         enddo
1481
1482         call chainbuild
1483         call overlap_sc_list(ioverlap,ioverlap_last)
1484 !        write (iout,*) 'Overlaping residues ',ioverlap_last,
1485 !     &           (ioverlap(j),j=1,ioverlap_last)
1486       enddo
1487
1488       if (k.le.1000.and.ioverlap_last.eq.0) then
1489         scfail=.false.
1490         if (had_overlaps) then
1491           write (iout,*) '#OVERLAPing all corrected after ',k,&
1492                ' random generation'
1493         endif
1494       else
1495         scfail=.true.
1496         write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
1497         write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
1498       endif
1499
1500       return
1501
1502  999  continue
1503       write (iout,'(a30,i5,a12,i4)') &
1504                      '#OVERLAP FAIL in gen_side after',maxsi,&
1505                      'iter for RES',i
1506       scfail=.true.
1507       return
1508       end subroutine overlap_sc
1509 !-----------------------------------------------------------------------------
1510       subroutine overlap_sc_list(ioverlap,ioverlap_last)
1511       use calc_data
1512 !      implicit real*8 (a-h,o-z)
1513 !      include 'DIMENSIONS'
1514 !      include 'COMMON.GEO'
1515 !      include 'COMMON.LOCAL'
1516 !      include 'COMMON.IOUNITS'
1517 !      include 'COMMON.CHAIN'
1518 !      include 'COMMON.INTERACT'
1519 !      include 'COMMON.FFIELD'
1520 !      include 'COMMON.VAR'
1521 !      include 'COMMON.CALC'
1522       logical :: fail
1523       integer,dimension(nres) :: ioverlap !(maxres)
1524       integer :: ioverlap_last
1525 !el local variables
1526       integer :: ind,iint
1527       real(kind=8) :: redfac,sig        !rrij,sigsq,
1528       integer :: itypi,itypj,itypi1
1529       real(kind=8) :: xi,yi,zi,sig0ij,rcomp,rrij,rij_shift
1530       data redfac /0.5D0/
1531
1532       ioverlap_last=0
1533 ! Check for SC-SC overlaps and mark residues
1534 !      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
1535       ind=0
1536       do i=iatsc_s,iatsc_e
1537         if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) cycle
1538         if (molnum(i).eq.5) print *,"WTF",i,iatsc_s,iatsc_e
1539         if (molnum(i).eq.5) cycle
1540         itypi=iabs(itype(i,molnum(i)))
1541         itypi1=iabs(itype(i+1,1))
1542         xi=c(1,nres+i)
1543         yi=c(2,nres+i)
1544         zi=c(3,nres+i)
1545         dxi=dc_norm(1,nres+i)
1546         dyi=dc_norm(2,nres+i)
1547         dzi=dc_norm(3,nres+i)
1548         print *,i,itypi,"sc_move"
1549         dsci_inv=dsc_inv(itypi)
1550 !
1551        do iint=1,nint_gr(i)
1552          do j=istart(i,iint),iend(i,iint)
1553          if (itype(j,molnum(j)).eq.ntyp1_molec(molnum(j))) cycle
1554             ind=ind+1
1555             itypj=iabs(itype(j,molnum(j)))
1556         print *,j,itypj,"sc_move"
1557
1558             dscj_inv=dsc_inv(itypj)
1559             sig0ij=sigma(itypi,itypj)
1560             chi1=chi(itypi,itypj)
1561             chi2=chi(itypj,itypi)
1562             chi12=chi1*chi2
1563             chip1=chip(itypi)
1564             chip2=chip(itypj)
1565             chip12=chip1*chip2
1566             alf1=alp(itypi)
1567             alf2=alp(itypj)
1568             alf12=0.5D0*(alf1+alf2)
1569           if (j.gt.i+1) then
1570            rcomp=sigmaii(itypi,itypj)
1571           else
1572            rcomp=sigma(itypi,itypj)
1573           endif
1574 !         print '(2(a3,2i3),a3,2f10.5)',
1575 !     &        ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
1576 !     &        ,rcomp
1577             xj=c(1,nres+j)-xi
1578             yj=c(2,nres+j)-yi
1579             zj=c(3,nres+j)-zi
1580             dxj=dc_norm(1,nres+j)
1581             dyj=dc_norm(2,nres+j)
1582             dzj=dc_norm(3,nres+j)
1583             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1584             rij=dsqrt(rrij)
1585             call sc_angular
1586             sigsq=1.0D0/sigsq
1587             sig=sig0ij*dsqrt(sigsq)
1588             rij_shift=1.0D0/rij-sig+sig0ij
1589
1590 !t          if ( 1.0/rij .lt. redfac*rcomp .or.
1591 !t     &       rij_shift.le.0.0D0 ) then
1592             if ( rij_shift.le.0.0D0 ) then
1593 !d           write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
1594 !d     &     'overlap SC-SC: i=',i,' j=',j,
1595 !d     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
1596 !d     &     rcomp,1.0/rij,rij_shift
1597           ioverlap_last=ioverlap_last+1
1598           ioverlap(ioverlap_last)=i
1599           do k=1,ioverlap_last-1
1600            if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
1601           enddo
1602           ioverlap_last=ioverlap_last+1
1603           ioverlap(ioverlap_last)=j
1604           do k=1,ioverlap_last-1
1605            if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
1606           enddo
1607          endif
1608         enddo
1609        enddo
1610       enddo
1611       return
1612       end subroutine overlap_sc_list
1613 #endif
1614 !-----------------------------------------------------------------------------
1615 ! energy_p_new_barrier.F
1616 !-----------------------------------------------------------------------------
1617       subroutine sc_angular
1618 ! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
1619 ! om12. Called by ebp, egb, and egbv.
1620       use calc_data
1621 !      implicit none
1622 !      include 'COMMON.CALC'
1623 !      include 'COMMON.IOUNITS'
1624       erij(1)=xj*rij
1625       erij(2)=yj*rij
1626       erij(3)=zj*rij
1627       om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1628       om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1629       om12=dxi*dxj+dyi*dyj+dzi*dzj
1630       chiom12=chi12*om12
1631 ! Calculate eps1(om12) and its derivative in om12
1632       faceps1=1.0D0-om12*chiom12
1633       faceps1_inv=1.0D0/faceps1
1634       eps1=dsqrt(faceps1_inv)
1635 ! Following variable is eps1*deps1/dom12
1636       eps1_om12=faceps1_inv*chiom12
1637 ! diagnostics only
1638 !      faceps1_inv=om12
1639 !      eps1=om12
1640 !      eps1_om12=1.0d0
1641 !      write (iout,*) "om12",om12," eps1",eps1
1642 ! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
1643 ! and om12.
1644       om1om2=om1*om2
1645       chiom1=chi1*om1
1646       chiom2=chi2*om2
1647       facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1648 !      print *,"TUT?",om1*chiom1,facsig,om1,om2,om12
1649       sigsq=1.0D0-facsig*faceps1_inv
1650       sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
1651       sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
1652       sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
1653 ! diagnostics only
1654 !      sigsq=1.0d0
1655 !      sigsq_om1=0.0d0
1656 !      sigsq_om2=0.0d0
1657 !      sigsq_om12=0.0d0
1658 !      write (iout,*) "chiom1",chiom1," chiom2",chiom2," chiom12",chiom12
1659 !      write (iout,*) "faceps1",faceps1," faceps1_inv",faceps1_inv,
1660 !     &    " eps1",eps1
1661 ! Calculate eps2 and its derivatives in om1, om2, and om12.
1662       chipom1=chip1*om1
1663       chipom2=chip2*om2
1664       chipom12=chip12*om12
1665       facp=1.0D0-om12*chipom12
1666       facp_inv=1.0D0/facp
1667       facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1668 !      write (iout,*) "chipom1",chipom1," chipom2",chipom2,
1669 !     &  " chipom12",chipom12," facp",facp," facp_inv",facp_inv
1670 ! Following variable is the square root of eps2
1671       eps2rt=1.0D0-facp1*facp_inv
1672 ! Following three variables are the derivatives of the square root of eps
1673 ! in om1, om2, and om12.
1674       eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
1675       eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
1676       eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
1677 ! Evaluate the "asymmetric" factor in the VDW constant, eps3
1678       eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
1679 !      write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
1680 !      write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
1681 !     &  " eps2rt_om12",eps2rt_om12
1682 ! Calculate whole angle-dependent part of epsilon and contributions
1683 ! to its derivatives
1684       return
1685       end subroutine sc_angular
1686 !-----------------------------------------------------------------------------
1687 ! initialize_p.F
1688       subroutine sc_angular_nucl
1689 ! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
1690 ! om12. Called by ebp, egb, and egbv.
1691 !      use calc_data
1692 !      implicit none
1693 !      include 'COMMON.CALC'
1694 !      include 'COMMON.IOUNITS'
1695       use comm_locel
1696       use calc_data_nucl
1697       erij(1)=xj*rij
1698       erij(2)=yj*rij
1699       erij(3)=zj*rij
1700       om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1701       om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1702       om12=dxi*dxj+dyi*dyj+dzi*dzj
1703       chiom12=chi12*om12
1704 ! Calculate eps1(om12) and its derivative in om12
1705       faceps1=1.0D0-om12*chiom12
1706       faceps1_inv=1.0D0/faceps1
1707       eps1=dsqrt(faceps1_inv)
1708 ! Following variable is eps1*deps1/dom12
1709       eps1_om12=faceps1_inv*chiom12
1710 ! diagnostics only
1711 !      faceps1_inv=om12
1712 !      eps1=om12
1713 !      eps1_om12=1.0d0
1714 !      write (iout,*) "om12",om12," eps1",eps1
1715 ! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
1716 ! and om12.
1717       om1om2=om1*om2
1718       chiom1=chi1*om1
1719       chiom2=chi2*om2
1720       facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1721       sigsq=1.0D0-facsig*faceps1_inv
1722       sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
1723       sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
1724       sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
1725       chipom1=chip1*om1
1726       chipom2=chip2*om2
1727       chipom12=chip12*om12
1728       facp=1.0D0-om12*chipom12
1729       facp_inv=1.0D0/facp
1730       facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1731 !      write (iout,*) "chipom1",chipom1," chipom2",chipom2,
1732 !     &  " chipom12",chipom12," facp",facp," facp_inv",facp_inv
1733 ! Following variable is the square root of eps2
1734       eps2rt=1.0D0-facp1*facp_inv
1735 ! Following three variables are the derivatives of the square root of eps
1736 ! in om1, om2, and om12.
1737       eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
1738       eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
1739       eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
1740 ! Evaluate the "asymmetric" factor in the VDW constant, eps3
1741       eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
1742 !      write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
1743 !      write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
1744 !     &  " eps2rt_om12",eps2rt_om12
1745 ! Calculate whole angle-dependent part of epsilon and contributions
1746 ! to its derivatives
1747       return
1748       end subroutine sc_angular_nucl
1749
1750 !-----------------------------------------------------------------------------
1751       subroutine int_bounds(total_ints,lower_bound,upper_bound)
1752 !      implicit real*8 (a-h,o-z)
1753 !      include 'DIMENSIONS'
1754       include 'mpif.h'
1755 !      include 'COMMON.SETUP'
1756       integer :: total_ints,lower_bound,upper_bound,nint
1757       integer,dimension(0:nfgtasks) :: int4proc,sint4proc       !(0:max_fg_procs)
1758       integer :: i,nexcess
1759       nint=total_ints/nfgtasks
1760       do i=1,nfgtasks
1761         int4proc(i-1)=nint
1762       enddo
1763       nexcess=total_ints-nint*nfgtasks
1764       do i=1,nexcess
1765         int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1766       enddo
1767       lower_bound=0
1768       do i=0,fg_rank-1
1769         lower_bound=lower_bound+int4proc(i)
1770       enddo
1771       upper_bound=lower_bound+int4proc(fg_rank)
1772       lower_bound=lower_bound+1
1773       return
1774       end subroutine int_bounds
1775 !-----------------------------------------------------------------------------
1776       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1777 !      implicit real*8 (a-h,o-z)
1778 !      include 'DIMENSIONS'
1779       include 'mpif.h'
1780 !      include 'COMMON.SETUP'
1781       integer :: total_ints,lower_bound,upper_bound,nint
1782       integer :: nexcess,i
1783       integer,dimension(0:nfgtasks) :: int4proc,sint4proc       !(0:max_fg_procs)
1784       nint=total_ints/nfgtasks1
1785       do i=1,nfgtasks1
1786         int4proc(i-1)=nint
1787       enddo
1788       nexcess=total_ints-nint*nfgtasks1
1789       do i=1,nexcess
1790         int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1791       enddo
1792       lower_bound=0
1793       do i=0,fg_rank1-1
1794         lower_bound=lower_bound+int4proc(i)
1795       enddo
1796       upper_bound=lower_bound+int4proc(fg_rank1)
1797       lower_bound=lower_bound+1
1798       return
1799       end subroutine int_bounds1
1800 !-----------------------------------------------------------------------------
1801 ! intcartderiv.F
1802 !-----------------------------------------------------------------------------
1803       subroutine chainbuild_cart
1804 !      implicit real*8 (a-h,o-z)
1805 !      include 'DIMENSIONS'
1806       use control_data
1807 #ifdef MPI
1808       include 'mpif.h'
1809 #endif
1810 !      include 'COMMON.SETUP'
1811 !      include 'COMMON.CHAIN'
1812 !      include 'COMMON.LOCAL'
1813 !      include 'COMMON.TIME1'
1814 !      include 'COMMON.IOUNITS'
1815       integer :: j,i,ierror,ierr
1816       real(kind=8) :: time00,time01
1817 #ifdef MPI
1818       if (nfgtasks.gt.1) then
1819 !        write (iout,*) "BCAST in chainbuild_cart"
1820 !        call flush(iout)
1821 ! Broadcast the order to build the chain and compute internal coordinates
1822 ! to the slaves. The slaves receive the order in ERGASTULUM.
1823         time00=MPI_Wtime()
1824 !      write (iout,*) "CHAINBUILD_CART: DC before BCAST"
1825 !      do i=0,nres
1826 !        write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1827 !     &   (dc(j,i+nres),j=1,3)
1828 !      enddo
1829         if (fg_rank.eq.0) &
1830           call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
1831         time_bcast7=time_bcast7+MPI_Wtime()-time00
1832         time01=MPI_Wtime()
1833         call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,&
1834           king,FG_COMM,IERR)
1835 !      write (iout,*) "CHAINBUILD_CART: DC after BCAST"
1836 !      do i=0,nres
1837 !        write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1838 !     &   (dc(j,i+nres),j=1,3)
1839 !      enddo
1840 !        write (iout,*) "End BCAST in chainbuild_cart"
1841 !        call flush(iout)
1842         time_bcast=time_bcast+MPI_Wtime()-time00
1843         time_bcastc=time_bcastc+MPI_Wtime()-time01
1844       endif
1845 #endif
1846       do j=1,3
1847         c(j,1)=dc(j,0)
1848       enddo
1849       do i=2,nres
1850         do j=1,3
1851           c(j,i)=c(j,i-1)+dc(j,i-1)
1852         enddo
1853       enddo
1854       do i=1,nres
1855         do j=1,3
1856           c(j,i+nres)=c(j,i)+dc(j,i+nres)
1857         enddo
1858       enddo
1859 !      write (iout,*) "CHAINBUILD_CART"
1860 !      call cartprint
1861       call int_from_cart1(.false.)
1862       return
1863       end subroutine chainbuild_cart
1864 !-----------------------------------------------------------------------------
1865 ! intcor.f
1866 !-----------------------------------------------------------------------------
1867       real(kind=8) function alpha(i1,i2,i3)
1868 !
1869 !  Calculates the planar angle between atoms (i1), (i2), and (i3).
1870 !
1871 !      implicit real*8 (a-h,o-z)
1872 !      include 'DIMENSIONS'
1873 !      include 'COMMON.GEO'
1874 !      include 'COMMON.CHAIN'
1875 !el local variables
1876       integer :: i1,i2,i3
1877       real(kind=8) :: x12,x23,y12,y23,z12,z23,vnorm,wnorm,scalar
1878       x12=c(1,i1)-c(1,i2)
1879       x23=c(1,i3)-c(1,i2)
1880       y12=c(2,i1)-c(2,i2)
1881       y23=c(2,i3)-c(2,i2)
1882       z12=c(3,i1)-c(3,i2)
1883       z23=c(3,i3)-c(3,i2)
1884       vnorm=dsqrt(x12*x12+y12*y12+z12*z12)
1885       wnorm=dsqrt(x23*x23+y23*y23+z23*z23)
1886       scalar=(x12*x23+y12*y23+z12*z23)/(vnorm*wnorm)
1887       alpha=arcos(scalar)
1888       return
1889       end function alpha
1890 !-----------------------------------------------------------------------------
1891       real(kind=8) function beta(i1,i2,i3,i4)
1892 !
1893 !  Calculates the dihedral angle between atoms (i1), (i2), (i3) and (i4)
1894 !
1895 !      implicit real*8 (a-h,o-z)
1896 !      include 'DIMENSIONS'
1897 !      include 'COMMON.GEO'
1898 !      include 'COMMON.CHAIN'
1899 !el local variables
1900       integer :: i1,i2,i3,i4
1901       real(kind=8) :: x12,x23,x34,y12,y23,y34,z12,z23,z34
1902       real(kind=8) :: wx,wy,wz,wnorm,vx,vy,vz,vnorm,scalar,angle
1903       real(kind=8) :: tx,ty,tz
1904       x12=c(1,i1)-c(1,i2)
1905       x23=c(1,i3)-c(1,i2)
1906       x34=c(1,i4)-c(1,i3)
1907       y12=c(2,i1)-c(2,i2)
1908       y23=c(2,i3)-c(2,i2)
1909       y34=c(2,i4)-c(2,i3)
1910       z12=c(3,i1)-c(3,i2)
1911       z23=c(3,i3)-c(3,i2)
1912       z34=c(3,i4)-c(3,i3)
1913 !d    print '(2i3,3f10.5)',i1,i2,x12,y12,z12
1914 !d    print '(2i3,3f10.5)',i2,i3,x23,y23,z23
1915 !d    print '(2i3,3f10.5)',i3,i4,x34,y34,z34
1916       wx=-y23*z34+y34*z23
1917       wy=x23*z34-z23*x34
1918       wz=-x23*y34+y23*x34
1919       wnorm=dsqrt(wx*wx+wy*wy+wz*wz)
1920       vx=y12*z23-z12*y23
1921       vy=-x12*z23+z12*x23
1922       vz=x12*y23-y12*x23
1923       vnorm=dsqrt(vx*vx+vy*vy+vz*vz)
1924       if (vnorm.gt.1.0D-13 .and. wnorm.gt.1.0D-13) then
1925       scalar=(vx*wx+vy*wy+vz*wz)/(vnorm*wnorm)
1926       if (dabs(scalar).gt.1.0D0) &
1927       scalar=0.99999999999999D0*scalar/dabs(scalar)
1928       angle=dacos(scalar)
1929 !d    print '(2i4,10f7.3)',i2,i3,vx,vy,vz,wx,wy,wz,vnorm,wnorm,
1930 !d   &scalar,angle
1931       else
1932       angle=pi
1933       endif
1934 !     if (angle.le.0.0D0) angle=pi+angle
1935       tx=vy*wz-vz*wy
1936       ty=-vx*wz+vz*wx
1937       tz=vx*wy-vy*wx
1938       scalar=tx*x23+ty*y23+tz*z23
1939       if (scalar.lt.0.0D0) angle=-angle
1940       beta=angle
1941       return
1942       end function beta
1943 !-----------------------------------------------------------------------------
1944       real(kind=8) function dist(i1,i2)
1945 !
1946 !  Calculates the distance between atoms (i1) and (i2).
1947 !
1948 !      implicit real*8 (a-h,o-z)
1949 !      include 'DIMENSIONS'
1950 !      include 'COMMON.GEO'
1951 !      include 'COMMON.CHAIN'
1952 !el local variables
1953       integer :: i1,i2
1954       real(kind=8) :: x12,y12,z12
1955       x12=c(1,i1)-c(1,i2)
1956       y12=c(2,i1)-c(2,i2)
1957       z12=c(3,i1)-c(3,i2)
1958       dist=dsqrt(x12*x12+y12*y12+z12*z12)
1959       return
1960       end function dist
1961 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
1962 !-----------------------------------------------------------------------------
1963 ! local_move.f
1964 !-----------------------------------------------------------------------------
1965       subroutine local_move_init(debug)
1966 !rc      implicit none
1967
1968 !     Includes
1969 !      implicit real*8 (a-h,o-z)
1970 !      include 'DIMENSIONS'  ! Needed by COMMON.LOCAL
1971 !      include 'COMMON.GEO'  ! For pi, deg2rad
1972 !      include 'COMMON.LOCAL'  ! For vbl
1973 !      include 'COMMON.LOCMOVE'
1974
1975 !     INPUT arguments
1976       logical :: debug
1977
1978
1979 !     Determine wheter to do some debugging output
1980       locmove_output=debug
1981
1982 !     Set the init_called flag to 1
1983       init_called=1
1984
1985 !     The following are never changed
1986       min_theta=60.D0*deg2rad  ! (0,PI)
1987       max_theta=175.D0*deg2rad  ! (0,PI)
1988       dmin2=vbl*vbl*2.*(1.-cos(min_theta))
1989       dmax2=vbl*vbl*2.*(1.-cos(max_theta))
1990       flag=1.0D300
1991       small=1.0D-5
1992       small2=0.5*small*small
1993
1994 !     Not really necessary...
1995       a_n=0
1996       b_n=0
1997       res_n=0
1998
1999       return
2000       end subroutine local_move_init
2001 !-----------------------------------------------------------------------------
2002       subroutine local_move(n_start, n_end, PHImin, PHImax)
2003 !     Perform a local move between residues m and n (inclusive)
2004 !     PHImin and PHImax [0,PI] determine the size of the move
2005 !     Works on whatever structure is in the variables theta and phi,
2006 !     sidechain variables are left untouched
2007 !     The final structure is NOT minimized, but both the cartesian
2008 !     variables c and the angles are up-to-date at the end (no further
2009 !     chainbuild is required)
2010 !rc      implicit none
2011       use random,only:ran_number
2012 !     Includes
2013 !      implicit real*8 (a-h,o-z)
2014 !      include 'DIMENSIONS'
2015 !      include 'COMMON.GEO'
2016 !      include 'COMMON.CHAIN'
2017 !      include 'COMMON.VAR'
2018 !      include 'COMMON.MINIM'
2019 !      include 'COMMON.SBRIDGE'
2020 !      include 'COMMON.LOCMOVE'
2021
2022 !     External functions
2023 !EL      integer move_res
2024 !EL      external move_res
2025 !EL      double precision ran_number
2026 !EL      external ran_number
2027
2028 !     INPUT arguments
2029       integer :: n_start, n_end  ! First and last residues to move
2030       real(kind=8) :: PHImin, PHImax  ! min/max angles [0,PI]
2031
2032 !     Local variables
2033       integer :: i,j
2034       real(kind=8) :: min,max
2035       integer :: iretcode
2036
2037
2038 !     Check if local_move_init was called.  This assumes that it
2039 !     would not be 1 if not explicitely initialized
2040       if (init_called.ne.1) then
2041         write(6,*)'   ***   local_move_init not called!!!'
2042         stop
2043       endif
2044
2045 !     Quick check for crazy range
2046       if (n_start.gt.n_end .or. n_start.lt.1 .or. n_end.gt.nres) then
2047         write(6,'(a,i3,a,i3)') &
2048              '   ***   Cannot make local move between n_start = ',&
2049              n_start,' and n_end = ',n_end
2050         return
2051       endif
2052
2053 !     Take care of end residues first...
2054       if (n_start.eq.1) then
2055 !     Move residue 1 (completely random)
2056         theta(3)=ran_number(min_theta,max_theta)
2057         phi(4)=ran_number(-PI,PI)
2058         i=2
2059       else
2060         i=n_start
2061       endif
2062       if (n_end.eq.nres) then
2063 !     Move residue nres (completely random)
2064         theta(nres)=ran_number(min_theta,max_theta)
2065         phi(nres)=ran_number(-PI,PI)
2066         j=nres-1
2067       else
2068         j=n_end
2069       endif
2070
2071 !     ...then go through all other residues one by one
2072 !     Start from the two extremes and converge
2073       call chainbuild
2074       do while (i.le.j)
2075         min=PHImin
2076         max=PHImax
2077 !$$$c     Move the first two residues by less than the others
2078 !$$$        if (i-n_start.lt.3) then
2079 !$$$          if (i-n_start.eq.0) then
2080 !$$$            min=0.4*PHImin
2081 !$$$            max=0.4*PHImax
2082 !$$$          else if (i-n_start.eq.1) then
2083 !$$$            min=0.8*PHImin
2084 !$$$            max=0.8*PHImax
2085 !$$$          else if (i-n_start.eq.2) then
2086 !$$$            min=PHImin
2087 !$$$            max=PHImax
2088 !$$$          endif
2089 !$$$        endif
2090
2091 !     The actual move, on residue i
2092         iretcode=move_res(min,max,i)  ! Discard iretcode
2093         i=i+1
2094
2095         if (i.le.j) then
2096           min=PHImin
2097           max=PHImax
2098 !$$$c     Move the last two residues by less than the others
2099 !$$$          if (n_end-j.lt.3) then
2100 !$$$            if (n_end-j.eq.0) then
2101 !$$$              min=0.4*PHImin
2102 !$$$              max=0.4*PHImax
2103 !$$$            else if (n_end-j.eq.1) then
2104 !$$$              min=0.8*PHImin
2105 !$$$              max=0.8*PHImax
2106 !$$$            else if (n_end-j.eq.2) then
2107 !$$$              min=PHImin
2108 !$$$              max=PHImax
2109 !$$$            endif
2110 !$$$          endif
2111
2112 !     The actual move, on residue j
2113           iretcode=move_res(min,max,j)  ! Discard iretcode
2114           j=j-1
2115         endif
2116       enddo
2117
2118       call int_from_cart(.false.,.false.)
2119
2120       return
2121       end subroutine local_move
2122 !-----------------------------------------------------------------------------
2123       subroutine output_tabs
2124 !     Prints out the contents of a_..., b_..., res_...
2125 !      implicit none
2126
2127 !     Includes
2128 !      include 'COMMON.GEO'
2129 !      include 'COMMON.LOCMOVE'
2130
2131 !     Local variables
2132       integer :: i,j
2133
2134       write(6,*)'a_...'
2135       write(6,'(8f7.1)')(a_ang(i)*rad2deg,i=0,a_n-1)
2136       write(6,'(8(2x,3l1,2x))')((a_tab(i,j),i=0,2),j=0,a_n-1)
2137
2138       write(6,*)'b_...'
2139       write(6,'(4f7.1)')(b_ang(i)*rad2deg,i=0,b_n-1)
2140       write(6,'(4(2x,3l1,2x))')((b_tab(i,j),i=0,2),j=0,b_n-1)
2141
2142       write(6,*)'res_...'
2143       write(6,'(12f7.1)')(res_ang(i)*rad2deg,i=0,res_n-1)
2144       write(6,'(12(2x,3l1,2x))')((res_tab(0,i,j),i=0,2),j=0,res_n-1)
2145       write(6,'(12(2x,3l1,2x))')((res_tab(1,i,j),i=0,2),j=0,res_n-1)
2146       write(6,'(12(2x,3l1,2x))')((res_tab(2,i,j),i=0,2),j=0,res_n-1)
2147
2148       return
2149       end subroutine output_tabs
2150 !-----------------------------------------------------------------------------
2151       subroutine angles2tab(PHImin,PHImax,n,ang,tab)
2152 !     Only uses angles if [0,PI] (but PHImin cannot be 0.,
2153 !     and PHImax cannot be PI)
2154 !      implicit none
2155
2156 !     Includes
2157 !      include 'COMMON.GEO'
2158
2159 !     INPUT arguments
2160       real(kind=8) :: PHImin,PHImax
2161
2162 !     OUTPUT arguments
2163       integer :: n
2164       real(kind=8),dimension(0:3) :: ang
2165       logical,dimension(0:2,0:3) :: tab
2166
2167
2168       if (PHImin .eq. PHImax) then
2169 !     Special case with two 010's
2170         n = 2;
2171         ang(0) = -PHImin;
2172         ang(1) = PHImin;
2173         tab(0,0) = .false.
2174         tab(2,0) = .false.
2175         tab(0,1) = .false.
2176         tab(2,1) = .false.
2177         tab(1,0) = .true.
2178         tab(1,1) = .true.
2179       else if (PHImin .eq. PI) then
2180 !     Special case with one 010
2181         n = 1
2182         ang(0) = PI
2183         tab(0,0) = .false.
2184         tab(2,0) = .false.
2185         tab(1,0) = .true.
2186       else if (PHImax .eq. 0.) then
2187 !     Special case with one 010
2188         n = 1
2189         ang(0) = 0.
2190         tab(0,0) = .false.
2191         tab(2,0) = .false.
2192         tab(1,0) = .true.
2193       else
2194 !     Standard cases
2195         n = 0
2196         if (PHImin .gt. 0.) then
2197 !     Start of range (011)
2198           ang(n) = PHImin
2199           tab(0,n) = .false.
2200           tab(1,n) = .true.
2201           tab(2,n) = .true.
2202 !     End of range (110)
2203           ang(n+1) = -PHImin
2204           tab(0,n+1) = .true.
2205           tab(1,n+1) = .true.
2206           tab(2,n+1) = .false.
2207           n = n+2
2208         endif
2209         if (PHImax .lt. PI) then
2210 !     Start of range (011)
2211           ang(n) = -PHImax
2212           tab(0,n) = .false.
2213           tab(1,n) = .true.
2214           tab(2,n) = .true.
2215 !     End of range (110)
2216           ang(n+1) = PHImax
2217           tab(0,n+1) = .true.
2218           tab(1,n+1) = .true.
2219           tab(2,n+1) = .false.
2220           n = n+2
2221         endif
2222       endif
2223
2224       return
2225       end subroutine angles2tab
2226 !-----------------------------------------------------------------------------
2227       subroutine minmax_angles(x,y,z,r,n,ang,tab)
2228 !     When solutions do not exist, assume all angles
2229 !     are acceptable - i.e., initial geometry must be correct
2230 !      implicit none
2231
2232 !     Includes
2233 !      include 'COMMON.GEO'
2234 !      include 'COMMON.LOCMOVE'
2235
2236 !     Input arguments
2237       real(kind=8) :: x,y,z,r
2238
2239 !     Output arguments
2240       integer :: n
2241       real(kind=8),dimension(0:3) :: ang
2242       logical,dimension(0:2,0:3) :: tab
2243
2244 !     Local variables
2245       real(kind=8) :: num, denom, phi
2246       real(kind=8) :: Kmin, Kmax
2247       integer :: i
2248
2249
2250       num = x*x + y*y + z*z
2251       denom = x*x + y*y
2252       n = 0
2253       if (denom .gt. 0.) then
2254         phi = atan2(y,x)
2255         denom = 2.*r*sqrt(denom)
2256         num = num+r*r
2257         Kmin = (num - dmin2)/denom
2258         Kmax = (num - dmax2)/denom
2259
2260 !     Allowed values of K (else all angles are acceptable)
2261 !     -1 <= Kmin <  1
2262 !     -1 <  Kmax <= 1
2263         if (Kmin .gt. 1. .or. abs(Kmin-1.) .lt. small2) then
2264           Kmin = -flag
2265         else if (Kmin .lt. -1. .or. abs(Kmin+1.) .lt. small2) then
2266           Kmin = PI
2267         else
2268           Kmin = acos(Kmin)
2269         endif
2270
2271         if (Kmax .lt. -1. .or. abs(Kmax+1.) .lt. small2) then
2272           Kmax = flag
2273         else if (Kmax .gt. 1. .or. abs(Kmax-1.) .lt. small2) then
2274           Kmax = 0.
2275         else
2276           Kmax = acos(Kmax)
2277         endif
2278
2279         if (Kmax .lt. Kmin) Kmax = Kmin
2280
2281         call angles2tab(Kmin, Kmax, n, ang, tab)
2282
2283 !     Add phi and check that angles are within range (-PI,PI]
2284         do i=0,n-1
2285           ang(i) = ang(i)+phi
2286           if (ang(i) .le. -PI) then
2287             ang(i) = ang(i)+2.*PI
2288           else if (ang(i) .gt. PI) then
2289             ang(i) = ang(i)-2.*PI
2290           endif
2291         enddo
2292       endif
2293
2294       return
2295       end subroutine minmax_angles
2296 !-----------------------------------------------------------------------------
2297       subroutine construct_tab
2298 !     Take a_... and b_... values and produces the results res_...
2299 !     x_ang are assumed to be all different (diff > small)
2300 !     x_tab(1,i) must be 1 for all i (i.e., all x_ang are acceptable)
2301 !      implicit none
2302
2303 !     Includes
2304 !      include 'COMMON.LOCMOVE'
2305
2306 !     Local variables
2307       integer :: n_max,i,j,index
2308       logical :: done
2309       real(kind=8) :: phi
2310
2311
2312       n_max = a_n + b_n
2313       if (n_max .eq. 0) then
2314         res_n = 0
2315         return
2316       endif
2317
2318       do i=0,n_max-1
2319         do j=0,1
2320           res_tab(j,0,i) = .true.
2321           res_tab(j,2,i) = .true.
2322           res_tab(j,1,i) = .false.
2323         enddo
2324       enddo
2325
2326       index = 0
2327       phi = -flag
2328       done = .false.
2329       do while (.not.done)
2330         res_ang(index) = flag
2331
2332 !     Check a first...
2333         do i=0,a_n-1
2334           if ((a_ang(i)-phi).gt.small .and. &
2335                a_ang(i) .lt. res_ang(index)) then
2336 !     Found a lower angle
2337             res_ang(index) = a_ang(i)
2338 !     Copy the values from a_tab into res_tab(0,,)
2339             res_tab(0,0,index) = a_tab(0,i)
2340             res_tab(0,1,index) = a_tab(1,i)
2341             res_tab(0,2,index) = a_tab(2,i)
2342 !     Set default values for res_tab(1,,)
2343             res_tab(1,0,index) = .true.
2344             res_tab(1,1,index) = .false.
2345             res_tab(1,2,index) = .true.
2346           else if (abs(a_ang(i)-res_ang(index)).lt.small) then
2347 !     Found an equal angle (can only be equal to a b_ang)
2348             res_tab(0,0,index) = a_tab(0,i)
2349             res_tab(0,1,index) = a_tab(1,i)
2350             res_tab(0,2,index) = a_tab(2,i)
2351           endif
2352         enddo
2353 !     ...then check b
2354         do i=0,b_n-1
2355           if ((b_ang(i)-phi).gt.small .and. &
2356                b_ang(i) .lt. res_ang(index)) then
2357 !     Found a lower angle
2358             res_ang(index) = b_ang(i)
2359 !     Copy the values from b_tab into res_tab(1,,)
2360             res_tab(1,0,index) = b_tab(0,i)
2361             res_tab(1,1,index) = b_tab(1,i)
2362             res_tab(1,2,index) = b_tab(2,i)
2363 !     Set default values for res_tab(0,,)
2364             res_tab(0,0,index) = .true.
2365             res_tab(0,1,index) = .false.
2366             res_tab(0,2,index) = .true.
2367           else if (abs(b_ang(i)-res_ang(index)).lt.small) then
2368 !     Found an equal angle (can only be equal to an a_ang)
2369             res_tab(1,0,index) = b_tab(0,i)
2370             res_tab(1,1,index) = b_tab(1,i)
2371             res_tab(1,2,index) = b_tab(2,i)
2372           endif
2373         enddo
2374
2375         if (res_ang(index) .eq. flag) then
2376           res_n = index
2377           done = .true.
2378         else if (index .eq. n_max-1) then
2379           res_n = n_max
2380           done = .true.
2381         else
2382           phi = res_ang(index)  ! Store previous angle
2383           index = index+1
2384         endif
2385       enddo
2386
2387 !     Fill the gaps
2388 !     First a...
2389       index = 0
2390       if (a_n .gt. 0) then
2391         do while (.not.res_tab(0,1,index))
2392           index=index+1
2393         enddo
2394         done = res_tab(0,2,index)
2395         do i=index+1,res_n-1
2396           if (res_tab(0,1,i)) then
2397             done = res_tab(0,2,i)
2398           else
2399             res_tab(0,0,i) = done
2400             res_tab(0,1,i) = done
2401             res_tab(0,2,i) = done
2402           endif
2403         enddo
2404         done = res_tab(0,0,index)
2405         do i=index-1,0,-1
2406           if (res_tab(0,1,i)) then
2407             done = res_tab(0,0,i)
2408           else
2409             res_tab(0,0,i) = done
2410             res_tab(0,1,i) = done
2411             res_tab(0,2,i) = done
2412           endif
2413         enddo
2414       else
2415         do i=0,res_n-1
2416           res_tab(0,0,i) = .true.
2417           res_tab(0,1,i) = .true.
2418           res_tab(0,2,i) = .true.
2419         enddo
2420       endif
2421 !     ...then b
2422       index = 0
2423       if (b_n .gt. 0) then
2424         do while (.not.res_tab(1,1,index))
2425           index=index+1
2426         enddo
2427         done = res_tab(1,2,index)
2428         do i=index+1,res_n-1
2429           if (res_tab(1,1,i)) then
2430             done = res_tab(1,2,i)
2431           else
2432             res_tab(1,0,i) = done
2433             res_tab(1,1,i) = done
2434             res_tab(1,2,i) = done
2435           endif
2436         enddo
2437         done = res_tab(1,0,index)
2438         do i=index-1,0,-1
2439           if (res_tab(1,1,i)) then
2440             done = res_tab(1,0,i)
2441           else
2442             res_tab(1,0,i) = done
2443             res_tab(1,1,i) = done
2444             res_tab(1,2,i) = done
2445           endif
2446         enddo
2447       else
2448         do i=0,res_n-1
2449           res_tab(1,0,i) = .true.
2450           res_tab(1,1,i) = .true.
2451           res_tab(1,2,i) = .true.
2452         enddo
2453       endif
2454
2455 !     Finally fill the last row with AND operation
2456       do i=0,res_n-1
2457         do j=0,2
2458           res_tab(2,j,i) = (res_tab(0,j,i) .and. res_tab(1,j,i))
2459         enddo
2460       enddo
2461
2462       return
2463       end subroutine construct_tab
2464 !-----------------------------------------------------------------------------
2465       subroutine construct_ranges(phi_n,phi_start,phi_end)
2466 !     Given the data in res_..., construct a table of
2467 !     min/max allowed angles
2468 !      implicit none
2469
2470 !     Includes
2471 !      include 'COMMON.GEO'
2472 !      include 'COMMON.LOCMOVE'
2473
2474 !     Output arguments
2475       integer :: phi_n
2476       real(kind=8),dimension(0:11) :: phi_start,phi_end
2477
2478 !     Local variables
2479       logical :: done
2480       integer :: index
2481
2482
2483       if (res_n .eq. 0) then
2484 !     Any move is allowed
2485         phi_n = 1
2486         phi_start(0) = -PI
2487         phi_end(0) = PI
2488       else
2489         phi_n = 0
2490         index = 0
2491         done = .false.
2492         do while (.not.done)
2493 !     Find start of range (01x)
2494           done = .false.
2495           do while (.not.done)
2496             if (res_tab(2,0,index).or.(.not.res_tab(2,1,index))) then
2497               index=index+1
2498             else
2499               done = .true.
2500               phi_start(phi_n) = res_ang(index)
2501             endif
2502             if (index .eq. res_n) done = .true.
2503           enddo
2504 !     If a start was found (index < res_n), find the end of range (x10)
2505 !     It may not be found without wrapping around
2506           if (index .lt. res_n) then
2507             done = .false.
2508             do while (.not.done)
2509               if ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) then
2510                 index=index+1
2511               else
2512                 done = .true.
2513               endif
2514               if (index .eq. res_n) done = .true.
2515             enddo
2516             if (index .lt. res_n) then
2517 !     Found the end of the range
2518               phi_end(phi_n) = res_ang(index)
2519               phi_n=phi_n+1
2520               index=index+1
2521               if (index .eq. res_n) then
2522                 done = .true.
2523               else
2524                 done = .false.
2525               endif
2526             else
2527 !     Need to wrap around
2528               done = .true.
2529               phi_end(phi_n) = flag
2530             endif
2531           endif
2532         enddo
2533 !     Take care of the last one if need to wrap around
2534         if (phi_end(phi_n) .eq. flag) then
2535           index = 0
2536           do while ((.not.res_tab(2,1,index)).or.res_tab(2,2,index))
2537             index=index+1
2538           enddo
2539           phi_end(phi_n) = res_ang(index) + 2.*PI
2540           phi_n=phi_n+1
2541         endif
2542       endif
2543
2544       return
2545       end subroutine construct_ranges
2546 !-----------------------------------------------------------------------------
2547       subroutine fix_no_moves(phi)
2548 !      implicit none
2549
2550 !     Includes
2551 !      include 'COMMON.GEO'
2552 !      include 'COMMON.LOCMOVE'
2553
2554 !     Output arguments
2555       real(kind=8) :: phi
2556
2557 !     Local variables
2558       integer :: index
2559       real(kind=8) :: diff,temp
2560
2561
2562 !     Look for first 01x in gammas (there MUST be at least one)
2563       diff = flag
2564       index = 0
2565       do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
2566         index=index+1
2567       enddo
2568       if (res_ang(index) .le. 0.D0) then ! Make sure it's from PHImax
2569 !     Try to increase PHImax
2570         if (index .gt. 0) then
2571           phi = res_ang(index-1)
2572           diff = abs(res_ang(index) - res_ang(index-1))
2573         endif
2574 !     Look for last (corresponding) x10
2575         index = res_n - 1
2576         do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
2577           index=index-1
2578         enddo
2579         if (index .lt. res_n-1) then
2580           temp = abs(res_ang(index) - res_ang(index+1))
2581           if (temp .lt. diff) then
2582             phi = res_ang(index+1)
2583             diff = temp
2584           endif
2585         endif
2586       endif
2587
2588 !     If increasing PHImax didn't work, decreasing PHImin
2589 !     will (with one exception)
2590 !     Look for first x10 (there MUST be at least one)
2591       index = 0
2592       do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
2593         index=index+1
2594       enddo
2595       if (res_ang(index) .lt. 0.D0) then ! Make sure it's from PHImin
2596 !     Try to decrease PHImin
2597         if (index .lt. res_n-1) then
2598           temp = abs(res_ang(index) - res_ang(index+1))
2599           if (res_ang(index+1) .le. 0.D0 .and. temp .lt. diff) then
2600             phi = res_ang(index+1)
2601             diff = temp
2602           endif
2603         endif
2604 !     Look for last (corresponding) 01x
2605         index = res_n - 1
2606         do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
2607           index=index-1
2608         enddo
2609         if (index .gt. 0) then
2610           temp = abs(res_ang(index) - res_ang(index-1))
2611           if (res_ang(index-1) .ge. 0.D0 .and. temp .lt. diff) then
2612             phi = res_ang(index-1)
2613             diff = temp
2614           endif
2615         endif
2616       endif
2617
2618 !     If it still didn't work, it must be PHImax == 0. or PHImin == PI
2619       if (diff .eq. flag) then
2620         index = 0
2621         if (res_tab(index,1,0) .or. (.not.res_tab(index,1,1)) .or. &
2622              res_tab(index,1,2)) index = res_n - 1
2623 !     This MUST work at this point
2624         if (index .eq. 0) then
2625           phi = res_ang(1)
2626         else
2627           phi = res_ang(index - 1)
2628         endif
2629       endif
2630
2631       return
2632       end subroutine fix_no_moves
2633 !-----------------------------------------------------------------------------
2634       integer function move_res(PHImin,PHImax,i_move)
2635 !     Moves residue i_move (in array c), leaving everything else fixed
2636 !     Starting geometry is not checked, it should be correct!
2637 !     R(,i_move) is the only residue that will move, but must have
2638 !     1 < i_move < nres (i.e., cannot move ends)
2639 !     Whether any output is done is controlled by locmove_output
2640 !rc      implicit none
2641       use random,only:ran_number
2642 !     Includes
2643 !      implicit real*8 (a-h,o-z)
2644 !      include 'DIMENSIONS'
2645 !      include 'COMMON.CHAIN'
2646 !      include 'COMMON.GEO'
2647 !      include 'COMMON.LOCMOVE'
2648
2649 !     External functions
2650 !EL      double precision ran_number
2651 !EL      external ran_number
2652
2653 !     Input arguments
2654       real(kind=8) :: PHImin,PHImax
2655       integer :: i_move
2656
2657 !     RETURN VALUES:
2658 !     0: move successfull
2659 !     1: Dmin or Dmax had to be modified
2660 !     2: move failed - check your input geometry
2661
2662
2663 !     Local variables
2664       real(kind=8),dimension(0:2) :: X,Y,Z,Orig
2665       real(kind=8),dimension(0:2) :: P
2666       logical :: no_moves,done
2667       integer :: index,i,j
2668       real(kind=8) :: phi,temp,radius
2669       real(kind=8),dimension(0:11) :: phi_start,phi_end
2670       integer :: phi_n
2671
2672 !     Set up the coordinate system
2673       do i=0,2
2674         Orig(i)=0.5*(c(i+1,i_move-1)+c(i+1,i_move+1)) ! Position of origin
2675       enddo
2676
2677       do i=0,2
2678         Z(i)=c(i+1,i_move+1)-c(i+1,i_move-1)
2679       enddo
2680       temp=sqrt(Z(0)*Z(0)+Z(1)*Z(1)+Z(2)*Z(2))
2681       do i=0,2
2682         Z(i)=Z(i)/temp
2683       enddo
2684
2685       do i=0,2
2686         X(i)=c(i+1,i_move)-Orig(i)
2687       enddo
2688 !     radius is the radius of the circle on which c(,i_move) can move
2689       radius=sqrt(X(0)*X(0)+X(1)*X(1)+X(2)*X(2))
2690       do i=0,2
2691         X(i)=X(i)/radius
2692       enddo
2693
2694       Y(0)=Z(1)*X(2)-X(1)*Z(2)
2695       Y(1)=X(0)*Z(2)-Z(0)*X(2)
2696       Y(2)=Z(0)*X(1)-X(0)*Z(1)
2697
2698 !     Calculate min, max angles coming from dmin, dmax to c(,i_move-2)
2699       if (i_move.gt.2) then
2700         do i=0,2
2701           P(i)=c(i+1,i_move-2)-Orig(i)
2702         enddo
2703         call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),&
2704              P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),&
2705              P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),&
2706              radius,a_n,a_ang,a_tab)
2707       else
2708         a_n=0
2709       endif
2710
2711 !     Calculate min, max angles coming from dmin, dmax to c(,i_move+2)
2712       if (i_move.lt.nres-2) then
2713         do i=0,2
2714           P(i)=c(i+1,i_move+2)-Orig(i)
2715         enddo
2716         call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),&
2717              P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),&
2718              P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),&
2719              radius,b_n,b_ang,b_tab)
2720       else
2721         b_n=0
2722       endif
2723
2724 !     Construct the resulting table for alpha and beta
2725       call construct_tab()
2726
2727       if (locmove_output) then
2728         print *,'ALPHAS & BETAS TABLE'
2729         call output_tabs()
2730       endif
2731
2732 !     Check that there is at least one possible move
2733       no_moves = .true.
2734       if (res_n .eq. 0) then
2735         no_moves = .false.
2736       else
2737         index = 0
2738         do while ((index .lt. res_n) .and. no_moves)
2739           if (res_tab(2,1,index)) no_moves = .false.
2740           index=index+1
2741         enddo
2742       endif
2743       if (no_moves) then
2744         if (locmove_output) print *,'   ***   Cannot move anywhere'
2745         move_res=2
2746         return
2747       endif
2748
2749 !     Transfer res_... into a_...
2750       a_n = 0
2751       do i=0,res_n-1
2752         if ( (res_tab(2,0,i).neqv.res_tab(2,1,i)) .or. &
2753              (res_tab(2,0,i).neqv.res_tab(2,2,i)) ) then
2754           a_ang(a_n) = res_ang(i)
2755           do j=0,2
2756             a_tab(j,a_n) = res_tab(2,j,i)
2757           enddo
2758           a_n=a_n+1
2759         endif
2760       enddo
2761
2762 !     Check that the PHI's are within [0,PI]
2763       if (PHImin .lt. 0. .or. abs(PHImin) .lt. small) PHImin = -flag
2764       if (PHImin .gt. PI .or. abs(PHImin-PI) .lt. small) PHImin = PI
2765       if (PHImax .gt. PI .or. abs(PHImax-PI) .lt. small) PHImax = flag
2766       if (PHImax .lt. 0. .or. abs(PHImax) .lt. small) PHImax = 0.
2767       if (PHImax .lt. PHImin) PHImax = PHImin
2768 !     Calculate min and max angles coming from PHImin and PHImax,
2769 !     and put them in b_...
2770       call angles2tab(PHImin, PHImax, b_n, b_ang, b_tab)
2771 !     Construct the final table
2772       call construct_tab()
2773
2774       if (locmove_output) then
2775         print *,'FINAL TABLE'
2776         call output_tabs()
2777       endif
2778
2779 !     Check that there is at least one possible move
2780       no_moves = .true.
2781       if (res_n .eq. 0) then
2782         no_moves = .false.
2783       else
2784         index = 0
2785         do while ((index .lt. res_n) .and. no_moves)
2786           if (res_tab(2,1,index)) no_moves = .false.
2787           index=index+1
2788         enddo
2789       endif
2790
2791       if (no_moves) then
2792 !     Take care of the case where no solution exists...
2793         call fix_no_moves(phi)
2794         if (locmove_output) then
2795           print *,'   ***   Had to modify PHImin or PHImax'
2796           print *,'phi: ',phi*rad2deg
2797         endif
2798         move_res=1
2799       else
2800 !     ...or calculate the solution
2801 !     Construct phi_start/phi_end arrays
2802         call construct_ranges(phi_n, phi_start, phi_end)
2803 !     Choose random angle phi in allowed range(s)
2804         temp = 0.
2805         do i=0,phi_n-1
2806           temp = temp + phi_end(i) - phi_start(i)
2807         enddo
2808         phi = ran_number(phi_start(0),phi_start(0)+temp)
2809         index = 0
2810         done = .false.
2811         do while (.not.done)
2812           if (phi .lt. phi_end(index)) then
2813             done = .true.
2814           else
2815             index=index+1
2816           endif
2817           if (index .eq. phi_n) then
2818             done = .true.
2819           else if (.not.done) then
2820             phi = phi + phi_start(index) - phi_end(index-1)
2821           endif
2822         enddo
2823         if (index.eq.phi_n) phi=phi_end(phi_n-1) ! Fix numerical errors
2824         if (phi .gt. PI) phi = phi-2.*PI
2825
2826         if (locmove_output) then
2827           print *,'ALLOWED RANGE(S)'
2828           do i=0,phi_n-1
2829             print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
2830           enddo
2831           print *,'phi: ',phi*rad2deg
2832         endif
2833         move_res=0
2834       endif
2835
2836 !     Re-use radius as temp variable
2837       temp=radius*cos(phi)
2838       radius=radius*sin(phi)
2839       do i=0,2
2840         c(i+1,i_move)=Orig(i)+temp*X(i)+radius*Y(i)
2841       enddo
2842
2843       return
2844       end function move_res
2845 !-----------------------------------------------------------------------------
2846       subroutine loc_test
2847 !rc      implicit none
2848
2849 !     Includes
2850 !      implicit real*8 (a-h,o-z)
2851 !      include 'DIMENSIONS'
2852 !      include 'COMMON.GEO'
2853 !      include 'COMMON.LOCAL'
2854 !      include 'COMMON.LOCMOVE'
2855
2856 !     External functions
2857 !EL      integer move_res
2858 !EL      external move_res
2859
2860 !     Local variables
2861       integer :: i,j,imov
2862       integer :: phi_n
2863       real(kind=8),dimension(0:11) :: phi_start,phi_end
2864       real(kind=8) :: phi
2865       real(kind=8),dimension(0:2,0:5) :: R
2866
2867       locmove_output=.true.
2868
2869 !      call angles2tab(30.*deg2rad,70.*deg2rad,a_n,a_ang,a_tab)
2870 !      call angles2tab(80.*deg2rad,130.*deg2rad,b_n,b_ang,b_tab)
2871 !      call minmax_angles(0.D0,3.8D0,0.D0,3.8D0,b_n,b_ang,b_tab)
2872 !      call construct_tab
2873 !      call output_tabs
2874
2875 !      call construct_ranges(phi_n,phi_start,phi_end)
2876 !      do i=0,phi_n-1
2877 !        print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
2878 !      enddo
2879
2880 !      call fix_no_moves(phi)
2881 !      print *,'NO MOVES FOUND, BEST PHI IS',phi*rad2deg
2882
2883       R(0,0)=0.D0
2884       R(1,0)=0.D0
2885       R(2,0)=0.D0
2886       R(0,1)=0.D0
2887       R(1,1)=-cos(28.D0*deg2rad)
2888       R(2,1)=-0.5D0-sin(28.D0*deg2rad)
2889       R(0,2)=0.D0
2890       R(1,2)=0.D0
2891       R(2,2)=-0.5D0
2892       R(0,3)=cos(30.D0*deg2rad)
2893       R(1,3)=0.D0
2894       R(2,3)=0.D0
2895       R(0,4)=0.D0
2896       R(1,4)=0.D0
2897       R(2,4)=0.5D0
2898       R(0,5)=0.D0
2899       R(1,5)=cos(26.D0*deg2rad)
2900       R(2,5)=0.5D0+sin(26.D0*deg2rad)
2901       do i=1,5
2902         do j=0,2
2903           R(j,i)=vbl*R(j,i)
2904         enddo
2905       enddo
2906 !      i=move_res(R(0,1),0.D0*deg2rad,180.D0*deg2rad)
2907       imov=2
2908       i=move_res(0.D0*deg2rad,180.D0*deg2rad,imov)
2909       print *,'RETURNED ',i
2910       print *,(R(i,3)/vbl,i=0,2)
2911
2912       return
2913       end subroutine loc_test
2914 #endif
2915 !-----------------------------------------------------------------------------
2916 ! matmult.f
2917 !-----------------------------------------------------------------------------
2918       subroutine MATMULT(A1,A2,A3)
2919 !      implicit real*8 (a-h,o-z)
2920 !      include 'DIMENSIONS'
2921 !el local variables
2922       integer :: i,j,k
2923       real(kind=8) :: A3IJ
2924
2925       real(kind=8),DIMENSION(3,3) :: A1,A2,A3
2926       real(kind=8),DIMENSION(3,3) :: AI3
2927       DO 1 I=1,3
2928         DO 2 J=1,3
2929           A3IJ=0.0
2930           DO 3 K=1,3
2931     3       A3IJ=A3IJ+A1(I,K)*A2(K,J)
2932           AI3(I,J)=A3IJ
2933     2   CONTINUE
2934     1 CONTINUE
2935       DO 4 I=1,3
2936       DO 4 J=1,3
2937     4   A3(I,J)=AI3(I,J)
2938       return
2939       end subroutine MATMULT
2940 !-----------------------------------------------------------------------------
2941 ! readpdb.F
2942 !-----------------------------------------------------------------------------
2943       subroutine int_from_cart(lside,lprn)
2944 !      implicit real*8 (a-h,o-z)
2945 !      include 'DIMENSIONS'
2946       use control_data,only:out1file
2947 #ifdef MPI
2948       include "mpif.h"
2949 #endif
2950 !      include 'COMMON.LOCAL'
2951 !      include 'COMMON.VAR'
2952 !      include 'COMMON.CHAIN'
2953 !      include 'COMMON.INTERACT'
2954 !      include 'COMMON.IOUNITS'
2955 !      include 'COMMON.GEO'
2956 !      include 'COMMON.NAMES'
2957 !      include 'COMMON.CONTROL'
2958 !      include 'COMMON.SETUP'
2959       character(len=3) :: seq,res
2960 !      character*5 atom
2961       character(len=80) :: card
2962       real(kind=8),dimension(3,20) :: sccor
2963       integer :: i,j,iti !el  rescode,
2964       logical :: lside,lprn
2965       real(kind=8) :: di,cosfac,sinfac
2966       integer :: nres2
2967       nres2=2*nres
2968
2969       if(me.eq.king.or..not.out1file)then
2970        if (lprn) then
2971         write (iout,'(/a)') &
2972         'Internal coordinates calculated from crystal structure.'
2973         if (lside) then
2974           write (iout,'(8a)') '  Res  ','       dvb','     Theta',&
2975        '     Gamma','    Dsc_id','       Dsc','     Alpha',&
2976        '     Beta '
2977         else
2978           write (iout,'(4a)') '  Res  ','       dvb','     Theta',&
2979        '     Gamma'
2980         endif
2981        endif
2982       endif
2983       do i=1,nres-1
2984 !       if (molnum(i).ne.1) cycle
2985 !in wham      do i=1,nres
2986         iti=itype(i,1)
2987         if (((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
2988        (iti.ne.ntyp1  .and. itype(i+1,1).ne.ntyp1)).and.molnum(i).eq.1) then
2989           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
2990 !test          stop
2991         endif
2992 !#ifndef WHAM_RUN
2993         vbld(i+1)=dist(i,i+1)
2994         vbld_inv(i+1)=1.0d0/vbld(i+1)
2995 !#endif
2996         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
2997         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
2998       enddo
2999 !el -----
3000 !#ifdef WHAM_RUN
3001 !      if (itype(1,1).eq.ntyp1) then
3002 !        do j=1,3
3003 !          c(j,1)=c(j,2)+(c(j,3)-c(j,4))
3004 !        enddo
3005 !      endif
3006 !      if (itype(nres,1).eq.ntyp1) then
3007 !        do j=1,3
3008 !          c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
3009 !        enddo
3010 !      endif
3011 !#endif
3012 !      if (unres_pdb) then
3013 !        if (itype(1,1).eq.21) then
3014 !          theta(3)=90.0d0*deg2rad
3015 !          phi(4)=180.0d0*deg2rad
3016 !          vbld(2)=3.8d0
3017 !          vbld_inv(2)=1.0d0/vbld(2)
3018 !        endif
3019 !        if (itype(nres,1).eq.21) then
3020 !          theta(nres)=90.0d0*deg2rad
3021 !          phi(nres)=180.0d0*deg2rad
3022 !          vbld(nres)=3.8d0
3023 !          vbld_inv(nres)=1.0d0/vbld(2)
3024 !        endif
3025 !      endif
3026       if (lside) then
3027         do i=2,nres-1
3028           do j=1,3
3029             c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) &
3030            +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
3031 ! in wham            c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)
3032           enddo
3033           iti=itype(i,1)
3034           di=dist(i,nres+i)
3035 !#ifndef WHAM_RUN
3036 ! 10/03/12 Adam: Correction for zero SC-SC bond length
3037
3038           if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) &
3039            di=dsc(itype(i,molnum(i)))
3040           vbld(i+nres)=di
3041           if (itype(i,1).ne.10) then
3042             vbld_inv(i+nres)=1.0d0/di
3043           else
3044             vbld_inv(i+nres)=0.0d0
3045           endif
3046 !#endif
3047           if (iti.ne.10) then
3048             alph(i)=alpha(nres+i,i,nres2+2)
3049             omeg(i)=beta(nres+i,i,nres2+2,i+1)
3050           endif
3051           if (iti.ne.0) then
3052           if(me.eq.king.or..not.out1file)then
3053            if (lprn) &
3054            write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
3055            rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),&
3056            rad2deg*alph(i),rad2deg*omeg(i)
3057           endif
3058           else
3059           if(me.eq.king.or..not.out1file)then
3060            if (lprn) &
3061            write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
3062            rad2deg*theta(i),rad2deg*phi(i),dsc(iti+1),vbld(nres+i),&
3063            rad2deg*alph(i),rad2deg*omeg(i)
3064           endif
3065           endif
3066         enddo
3067       else if (lprn) then
3068         do i=2,nres
3069           iti=itype(i,1)
3070           if(me.eq.king.or..not.out1file) &
3071            write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,dist(i,i-1),&
3072            rad2deg*theta(i),rad2deg*phi(i)
3073         enddo
3074       endif
3075       return
3076       end subroutine int_from_cart
3077 !-----------------------------------------------------------------------------
3078       subroutine sc_loc_geom(lprn)
3079 !      implicit real*8 (a-h,o-z)
3080 !      include 'DIMENSIONS'
3081       use control_data,only:out1file
3082 #ifdef MPI
3083       include "mpif.h"
3084 #endif
3085 !      include 'COMMON.LOCAL'
3086 !      include 'COMMON.VAR'
3087 !      include 'COMMON.CHAIN'
3088 !      include 'COMMON.INTERACT'
3089 !      include 'COMMON.IOUNITS'
3090 !      include 'COMMON.GEO'
3091 !      include 'COMMON.NAMES'
3092 !      include 'COMMON.CONTROL'
3093 !      include 'COMMON.SETUP'
3094       real(kind=8),dimension(3) :: x_prime,y_prime,z_prime
3095       logical :: lprn
3096 !el local variables
3097       integer :: i,j,it,iti
3098       real(kind=8) :: cosfac2,sinfac2,xx,yy,zz,cosfac,sinfac
3099       do i=1,nres-1
3100         do j=1,3
3101           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
3102         enddo
3103       enddo
3104       do i=2,nres-1
3105         if (itype(i,1).ne.10) then
3106           do j=1,3
3107             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
3108           enddo
3109         else
3110           do j=1,3
3111             dc_norm(j,i+nres)=0.0d0
3112           enddo
3113         endif
3114       enddo
3115       do i=2,nres-1
3116         costtab(i+1) =dcos(theta(i+1))
3117         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
3118         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
3119         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
3120         cosfac2=0.5d0/(1.0d0+costtab(i+1))
3121         cosfac=dsqrt(cosfac2)
3122         sinfac2=0.5d0/(1.0d0-costtab(i+1))
3123         sinfac=dsqrt(sinfac2)
3124         it=itype(i,1)
3125
3126         if ((it.ne.10).and.(it.ne.ntyp1)) then
3127 !el        if (it.ne.10) then
3128 !
3129 !  Compute the axes of tghe local cartesian coordinates system; store in
3130 !   x_prime, y_prime and z_prime
3131 !
3132         do j=1,3
3133           x_prime(j) = 0.00
3134           y_prime(j) = 0.00
3135           z_prime(j) = 0.00
3136         enddo
3137         do j = 1,3
3138           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
3139           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
3140         enddo
3141         call vecpr(x_prime,y_prime,z_prime)
3142 !
3143 ! Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
3144 ! to local coordinate system. Store in xx, yy, zz.
3145 !
3146         xx=0.0d0
3147         yy=0.0d0
3148         zz=0.0d0
3149         do j = 1,3
3150           xx = xx + x_prime(j)*dc_norm(j,i+nres)
3151           yy = yy + y_prime(j)*dc_norm(j,i+nres)
3152           zz = zz + z_prime(j)*dc_norm(j,i+nres)
3153         enddo
3154
3155         xxref(i)=xx
3156         yyref(i)=yy
3157         zzref(i)=zz
3158         else
3159         xxref(i)=0.0d0
3160         yyref(i)=0.0d0
3161         zzref(i)=0.0d0
3162         endif
3163       enddo
3164       if (lprn) then
3165         do i=2,nres
3166           iti=itype(i,1)
3167           if(me.eq.king.or..not.out1file) &
3168            write (iout,'(a3,i4,3f10.5)') restyp(iti,1),i,xxref(i),&
3169             yyref(i),zzref(i)
3170         enddo
3171       endif
3172
3173       return
3174       end subroutine sc_loc_geom
3175 !-----------------------------------------------------------------------------
3176       subroutine sccenter(ires,nscat,sccor)
3177 !      implicit real*8 (a-h,o-z)
3178 !      include 'DIMENSIONS'
3179 !      include 'COMMON.CHAIN'
3180       integer :: i,j,ires,nscat
3181       real(kind=8),dimension(3,20) :: sccor
3182       real(kind=8) :: sccmj
3183 !        print *,"I am in sccenter",ires,nscat
3184       do j=1,3
3185         sccmj=0.0D0
3186         do i=1,nscat
3187           sccmj=sccmj+sccor(j,i)
3188 !C          print *,"insccent", ires,sccor(j,i)
3189         enddo
3190         dc(j,ires)=sccmj/nscat
3191       enddo
3192       return
3193       end subroutine sccenter
3194 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3195 !-----------------------------------------------------------------------------
3196       subroutine bond_regular
3197       use calc_data
3198 !      implicit real*8 (a-h,o-z)
3199 !      include 'DIMENSIONS'
3200 !      include 'COMMON.VAR'
3201 !      include 'COMMON.LOCAL'
3202 !      include 'COMMON.CALC'
3203 !      include 'COMMON.INTERACT'
3204 !      include 'COMMON.CHAIN'
3205       do i=1,nres-1
3206
3207        vbld(i+1)=vbl
3208        vbld_inv(i+1)=1.0d0/vbld(i+1)
3209        vbld(i+1+nres)=dsc(itype(i+1,molnum(i)))
3210        vbld_inv(i+1+nres)=dsc_inv(itype(i+1,molnum(i)))
3211 !       print *,vbld(i+1),vbld(i+1+nres)
3212       enddo
3213       return
3214       end subroutine bond_regular
3215 #endif
3216 !-----------------------------------------------------------------------------
3217 ! refsys.f
3218 !-----------------------------------------------------------------------------
3219       subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
3220 ! This subroutine calculates unit vectors of a local reference system
3221 ! defined by atoms (i2), (i3), and (i4). The x axis is the axis from
3222 ! atom (i3) to atom (i2), and the xy plane is the plane defined by atoms
3223 ! (i2), (i3), and (i4). z axis is directed according to the sign of the
3224 ! vector product (i3)-(i2) and (i3)-(i4). Sets fail to .true. if atoms
3225 ! (i2) and (i3) or (i3) and (i4) coincide or atoms (i2), (i3), and (i4)
3226 ! form a linear fragment. Returns vectors e1, e2, and e3.
3227 !      implicit real*8 (a-h,o-z)
3228 !      include 'DIMENSIONS'
3229       logical :: fail
3230       real(kind=8),dimension(3) :: e1,e2,e3
3231       real(kind=8),dimension(3) :: u,z
3232 !      include 'COMMON.IOUNITS'
3233 !      include 'COMMON.CHAIN'
3234       real(kind=8) :: coinc=1.0D-13,align=1.0D-13
3235 !el local variables
3236       integer :: i,i1,i2,i3,i4
3237       real(kind=8) :: v1,v2,v3,s1,s2,zi,ui,anorm
3238       fail=.false.
3239       s1=0.0
3240       s2=0.0
3241       do 1 i=1,3
3242       zi=c(i,i2)-c(i,i3)
3243       ui=c(i,i4)-c(i,i3)
3244       s1=s1+zi*zi
3245       s2=s2+ui*ui
3246       z(i)=zi
3247     1 u(i)=ui
3248       s1=sqrt(s1)
3249       s2=sqrt(s2)
3250       if (s1.gt.coinc) goto 2
3251       write (iout,1000) i2,i3,i1
3252       fail=.true.
3253 !     do 3 i=1,3
3254 !   3 c(i,i1)=0.0D0
3255       return
3256     2 if (s2.gt.coinc) goto 4
3257       write(iout,1000) i3,i4,i1
3258       fail=.true.
3259       do 5 i=1,3
3260     5 c(i,i1)=0.0D0
3261       return
3262     4 s1=1.0/s1
3263       s2=1.0/s2
3264       v1=z(2)*u(3)-z(3)*u(2)
3265       v2=z(3)*u(1)-z(1)*u(3)
3266       v3=z(1)*u(2)-z(2)*u(1)
3267       anorm=dsqrt(v1*v1+v2*v2+v3*v3)
3268       if (anorm.gt.align) goto 6
3269       write (iout,1010) i2,i3,i4,i1
3270       fail=.true.
3271 !     do 7 i=1,3
3272 !   7 c(i,i1)=0.0D0
3273       return
3274     6 anorm=1.0D0/anorm
3275       e3(1)=v1*anorm
3276       e3(2)=v2*anorm
3277       e3(3)=v3*anorm
3278       e1(1)=z(1)*s1
3279       e1(2)=z(2)*s1
3280       e1(3)=z(3)*s1
3281       e2(1)=e1(3)*e3(2)-e1(2)*e3(3)
3282       e2(2)=e1(1)*e3(3)-e1(3)*e3(1)
3283       e2(3)=e1(2)*e3(1)-e1(1)*e3(2)
3284  1000 format (/1x,' * * * Error - atoms',i4,' and',i4,' coincide.',&
3285        'coordinates of atom',i4,' are set to zero.')
3286  1010 format (/1x,' * * * Error - atoms',2(i4,2h, ),i4,' form a linear',&
3287        ' fragment. coordinates of atom',i4,' are set to zero.')
3288       return
3289       end subroutine refsys
3290 !-----------------------------------------------------------------------------
3291 ! int_to_cart.f
3292 !-----------------------------------------------------------------------------
3293       subroutine int_to_cart
3294 !--------------------------------------------------------------
3295 !  This subroutine converts the energy derivatives from internal
3296 !  coordinates to cartesian coordinates
3297 !-------------------------------------------------------------
3298 !      implicit real*8 (a-h,o-z)
3299 !      include 'DIMENSIONS'
3300 !      include 'COMMON.VAR'
3301 !      include 'COMMON.CHAIN'
3302 !      include 'COMMON.DERIV'
3303 !      include 'COMMON.GEO'
3304 !      include 'COMMON.LOCAL'
3305 !      include 'COMMON.INTERACT'
3306 !      include 'COMMON.MD'
3307 !      include 'COMMON.IOUNITS'
3308 !      include 'COMMON.SCCOR'
3309 !   calculating dE/ddc1
3310 !el local variables
3311        integer :: j,i
3312 !       print *,"gloc",gloc(:,:)
3313 !       print *, "gcart",gcart(:,:)
3314        if (nres.lt.3) go to 18
3315        do j=1,3
3316          gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
3317            +gloc(nres-2,icg)*dtheta(j,1,3)
3318 !         write(iout,*) "pierwszy gcart", gcart(j,2)
3319           if ((itype(2,1).ne.10).and.&
3320           (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).ne.5))) then
3321           gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ &
3322           gloc(ialph(2,1)+nside,icg)*domega(j,1,2)
3323         endif
3324        enddo
3325 !     Calculating the remainder of dE/ddc2
3326        do j=1,3
3327          gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ &
3328          gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
3329 !         write(iout,*) "drugi gcart", gcart(j,2)
3330         if((itype(2,1).ne.10).and.(molnum(2).ne.5)) then
3331           gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ &
3332           gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
3333         endif
3334         if(itype(3,1).ne.10) then
3335           gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+ &
3336           gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
3337         endif
3338         if(nres.gt.4) then
3339           gcart(j,2)=gcart(j,2)+gloc(2,icg)*dphi(j,1,5)
3340         endif
3341        enddo
3342 !  If there are only five residues
3343        if(nres.eq.5) then
3344          do j=1,3
3345            gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* &
3346            dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* &
3347            dtheta(j,1,5)
3348 !         if(itype(3,1).ne.10) then
3349           if ((itype(3,1).ne.10).and.&
3350           (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) then
3351            gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* &
3352            dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
3353          endif
3354 !        if(itype(4,1).ne.10) then
3355           if ((itype(4,1).ne.10).and.&
3356           (itype(4,molnum(4)).ne.ntyp1_molec(molnum(4)))) then
3357            gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* &
3358            dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
3359          endif
3360         enddo
3361        endif
3362 !    If there are more than five residues
3363       if(nres.gt.5) then
3364         do i=3,nres-3
3365          do j=1,3
3366           gcart(j,i)=gcart(j,i)+gloc(i-2,icg)*dphi(j,3,i+1) &
3367           +gloc(i-1,icg)*dphi(j,2,i+2)+ &
3368           gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+ &
3369           gloc(nres+i-3,icg)*dtheta(j,1,i+2)
3370           if((itype(i,1).ne.10).and.(molnum(nres-1).ne.5)) then
3371            gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+ &
3372            gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
3373           endif
3374           if(itype(i+1,1).ne.10) then
3375            gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1) &
3376            +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
3377           endif
3378          enddo
3379         enddo
3380       endif
3381 !  Setting dE/ddnres-2
3382       if(nres.gt.5) then
3383          do j=1,3
3384            gcart(j,nres-2)=gcart(j,nres-2)+gloc(nres-4,icg)* &
3385            dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres) &
3386            +gloc(2*nres-6,icg)* &
3387            dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
3388           if((itype(nres-2,1).ne.10).and.(molnum(nres-1).ne.5)) then
3389               gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)* &
3390               dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)* &
3391               domega(j,2,nres-2)
3392           endif
3393           if((itype(nres-1,1).ne.10).and.(molnum(nres-1).ne.5)) then
3394              gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)* &
3395              dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
3396              domega(j,1,nres-1)
3397           endif
3398          enddo
3399       endif
3400 !  Settind dE/ddnres-1
3401 #define DEBUG
3402 #ifdef DEBUG
3403           j=1
3404               write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
3405         gloc(2*nres-5,icg),dtheta(j,2,nres)
3406
3407 #endif
3408 !#undef DEBUG
3409
3410        do j=1,3
3411         gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
3412         gloc(2*nres-5,icg)*dtheta(j,2,nres)
3413 !#define DEBUG
3414 #ifdef DEBUG
3415               write(iout,*)"in int to cartb",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
3416         gloc(2*nres-5,icg),dtheta(j,2,nres)
3417
3418 #endif
3419 !#undef DEBUG
3420         if((itype(nres-1,1).ne.10).and.(molnum(nres-1).ne.5)) then
3421           gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
3422           dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
3423           domega(j,2,nres-1)
3424 !#define DEBUG
3425 #ifdef DEBUG
3426               write(iout,*)"in int to cart2",i,gcart(j,nres-1),gloc(ialph(nres-1,1),icg)* &
3427           dalpha(j,2,nres-1),gloc(ialph(nres-1,1)+nside,icg), &
3428           domega(j,2,nres-1)
3429
3430 #endif
3431 !#undef DEBUG
3432
3433         endif
3434         enddo
3435 !   The side-chain vector derivatives
3436         do i=2,nres-1
3437          if(itype(i,1).ne.10 .and.  &
3438            itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)).and.&
3439              (molnum(i).ne.5)) then
3440             do j=1,3
3441               gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
3442               +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
3443 !#define DEBUG
3444 #ifdef DEBUG
3445               write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), &
3446               gloc(ialph(i,1)+nside,icg),domega(j,3,i)
3447 #endif
3448 #undef DEBUG
3449             enddo
3450          endif
3451        enddo
3452 !----------------------------------------------------------------------
3453 ! INTERTYP=1 SC...Ca...Ca...Ca
3454 ! INTERTYP=2 Ca...Ca...Ca...SC
3455 ! INTERTYP=3 SC...Ca...Ca...SC
3456 !   calculating dE/ddc1
3457   18   continue
3458 !         write(iout,*) "przed sccor gcart", gcart(1,2)
3459
3460 !       do i=1,nres
3461 !       gloc(i,icg)=0.0D0
3462 !          write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
3463 !       enddo
3464        if (nres.lt.2) return
3465        if ((nres.lt.3).and.(itype(1,1).eq.10)) return
3466        if ((itype(1,1).ne.10).and. &
3467         (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) then
3468         do j=1,3
3469 !c Derviative was calculated for oposite vector of side chain therefore
3470 ! there is "-" sign before gloc_sc
3471          gxcart(j,1)=gxcart(j,1)-gloc_sc(1,0,icg)* &
3472            dtauangle(j,1,1,3)
3473          gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
3474            dtauangle(j,1,2,3)
3475           if ((itype(2,1).ne.10).and. &
3476         (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2))).and.(molnum(2).ne.5)) then
3477          gxcart(j,1)= gxcart(j,1) &
3478                      -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
3479          gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
3480             dtauangle(j,3,2,3)
3481           endif
3482        enddo
3483        endif
3484          if ((nres.ge.3).and.(itype(3,molnum(3)).ne.10).and.&
3485          (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) &
3486       then
3487          do j=1,3
3488          gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
3489          enddo
3490          endif
3491 !   As potetnial DO NOT depend on omicron anlge their derivative is
3492 !   ommited
3493 !     &     +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3)
3494
3495 !         write(iout,*) "przed dE/ddc2 gcart", gcart(1,2)
3496
3497 !     Calculating the remainder of dE/ddc2
3498        do j=1,3
3499          if((itype(2,1).ne.10).and. &
3500            (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).ne.5))) then
3501            if ((itype(1,1).ne.10).and.&
3502               ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))).and.(molnum(1).ne.5))&
3503             gxcart(j,2)=gxcart(j,2)+ &
3504                                gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
3505         if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3)).and.molnum(3).ne.5) &
3506          then
3507            gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
3508 !c                  the   - above is due to different vector direction
3509            gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4)
3510           endif
3511           if (nres.gt.3) then
3512 !           if ((itype(1,1).ne.10).and.&
3513 !              ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))))) &
3514            gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4)
3515 !c                  the   - above is due to different vector direction
3516            gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4)
3517 !          write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart",gcart(j,2)
3518 !           write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
3519           endif
3520          endif
3521          if ((itype(1,1).ne.10).and.&
3522          (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) then
3523           gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
3524 !           write(iout,*)  gloc_sc(1,0,icg),dtauangle(j,1,3,3)
3525         endif
3526          if ((itype(3,1).ne.10).and.(nres.ge.3).and.(molnum(3).ne.5)) then
3527           gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
3528 !           write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
3529          endif
3530          if ((itype(4,1).ne.10).and.(nres.ge.4).and.(molnum(4).ne.5)) then
3531           gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
3532 !           write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
3533          endif
3534
3535 !      write(iout,*) gcart(j,2),itype(2,1),itype(1,1),itype(3,1), "gcart2"
3536        enddo
3537 !         write(iout,*) "po dE/ddc2 gcart", gcart(1,2)
3538
3539 !    If there are more than five residues
3540       if(nres.ge.5) then
3541         do i=3,nres-2
3542          do j=1,3
3543 !          write(iout,*) "before", gcart(j,i)
3544           if ((itype(i,1).ne.10).and.&
3545           (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))).and.(molnum(i).ne.5)) then
3546           gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
3547           *dtauangle(j,2,3,i+1) &
3548           -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
3549           gcart(j,i)=gcart(j,i)+gloc_sc(1,i-1,icg) &
3550           *dtauangle(j,1,2,i+2)
3551 !                   write(iout,*) "new",j,i,
3552 !     &  gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
3553 !          if (itype(i-1,1).ne.10) then
3554           if ((itype(i-1,1).ne.10).and.&
3555           (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1))).and.(molnum(i-1).eq.5)) then
3556
3557            gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
3558       *dtauangle(j,3,3,i+1)
3559           endif
3560 !          if (itype(i+1,1).ne.10) then
3561           if ((itype(i+1,1).ne.10).and.&
3562           (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1))).and.&
3563            (molnum(i+1).ne.5)) then
3564           gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
3565       *dtauangle(j,3,1,i+2)
3566            gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) &
3567       *dtauangle(j,3,2,i+2)
3568           endif
3569           endif
3570 !          if (itype(i-1,1).ne.10) then
3571           if ((itype(i-1,1).ne.10).and.&
3572           (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1))).and.&
3573          (molnum(i-1).ne.5)) then
3574            gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
3575            dtauangle(j,1,3,i+1)
3576           endif
3577 !          if (itype(i+1,1).ne.10) then
3578           if ((itype(i+1,1).ne.10).and.&
3579           (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))&
3580        .and. (molnum(i+1).ne.5)) then
3581            gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* &
3582            dtauangle(j,2,2,i+2)
3583 !          write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
3584 !     &    dtauangle(j,2,2,i+2)
3585           endif
3586 !          if (itype(i+2,1).ne.10) then
3587           if ((itype(i+2,1).ne.10).and.&
3588           (itype(i+2,molnum(i+2)).ne.ntyp1_molec(molnum(i+2))).and.(molnum(i+2).ne.5)) then
3589            gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
3590            dtauangle(j,2,1,i+3)
3591           endif
3592          enddo
3593         enddo
3594       endif
3595 !  Setting dE/ddnres-1
3596       if(nres.ge.4) then
3597          do j=1,3
3598          if ((itype(nres-1,1).ne.10).and.&
3599        (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1))).and.(molnum(nres-1).ne.5)) then
3600          gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
3601           *dtauangle(j,2,3,nres)
3602 !          write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
3603 !     &     dtauangle(j,2,3,nres), gxcart(j,nres-1)
3604 !         if (itype(nres-2,1).ne.10) then
3605          if ((itype(nres-2,1).ne.10).and.&
3606        (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).ne.5)) then
3607        gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
3608           *dtauangle(j,3,3,nres)
3609           endif
3610          if ((itype(nres,1).ne.10).and.&
3611          (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5)) then
3612         gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
3613           *dtauangle(j,3,1,nres+1)
3614         gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
3615           *dtauangle(j,3,2,nres+1)
3616           endif
3617          endif
3618          if ((itype(nres-2,1).ne.10).and.&
3619          (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).ne.5)) then
3620             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
3621          dtauangle(j,1,3,nres)
3622          endif
3623           if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5)) then
3624             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
3625            dtauangle(j,2,2,nres+1)
3626 !           write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
3627 !     &     dtauangle(j,2,2,nres+1), itype(nres-1,1),itype(nres,1)
3628            endif
3629          enddo
3630       endif
3631 !  Settind dE/ddnres
3632        if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
3633           (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5))then
3634        do j=1,3
3635         gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
3636        *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
3637        *dtauangle(j,2,3,nres+1)
3638         enddo
3639        endif
3640 !       write(iout,*) "final gcart",gcart(1,2)
3641 !   The side-chain vector derivatives
3642 !       print *,"gcart",gcart(:,:)
3643       return
3644       end subroutine int_to_cart
3645 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3646 !-----------------------------------------------------------------------------
3647 ! readrtns_CSA.F
3648 !-----------------------------------------------------------------------------
3649       subroutine gen_dist_constr
3650 ! Generate CA distance constraints.
3651 !      implicit real*8 (a-h,o-z)
3652 !      include 'DIMENSIONS'
3653 !      include 'COMMON.IOUNITS'
3654 !      include 'COMMON.GEO'
3655 !      include 'COMMON.VAR'
3656 !      include 'COMMON.INTERACT'
3657 !      include 'COMMON.LOCAL'
3658 !      include 'COMMON.NAMES'
3659 !      include 'COMMON.CHAIN'
3660 !      include 'COMMON.FFIELD'
3661 !      include 'COMMON.SBRIDGE'
3662 !      include 'COMMON.HEADER'
3663 !      include 'COMMON.CONTROL'
3664 !      include 'COMMON.DBASE'
3665 !      include 'COMMON.THREAD'
3666 !      include 'COMMON.TIME1'
3667 !      integer :: itype_pdb !(maxres)
3668 !      common /pizda/ itype_pdb(nres)
3669       character(len=2) :: iden
3670 !el local variables
3671       integer :: i,j
3672 !d      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
3673 !d      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
3674 !d     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
3675 !d     & ' nsup',nsup
3676       do i=nstart_sup,nstart_sup+nsup-1
3677 !d      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
3678 !d     &    ' seq_pdb', restyp(itype_pdb(i))
3679         do j=i+2,nstart_sup+nsup-1
3680           nhpb=nhpb+1
3681           ihpb(nhpb)=i+nstart_seq-nstart_sup
3682           jhpb(nhpb)=j+nstart_seq-nstart_sup
3683           forcon(nhpb)=weidis
3684           dhpb(nhpb)=dist(i,j)
3685         enddo
3686       enddo
3687 !d      write (iout,'(a)') 'Distance constraints:'
3688 !d      do i=nss+1,nhpb
3689 !d        ii=ihpb(i)
3690 !d        jj=jhpb(i)
3691 !d        iden='CA'
3692 !d        if (ii.gt.nres) then
3693 !d          iden='SC'
3694 !d          ii=ii-nres
3695 !d          jj=jj-nres
3696 !d        endif
3697 !d        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
3698 !d     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
3699 !d     &  dhpb(i),forcon(i)
3700 !d      enddo
3701 !      deallocate(itype_pdb)
3702
3703       return
3704       end subroutine gen_dist_constr
3705 #endif
3706 !-----------------------------------------------------------------------------
3707 ! cartprint.f
3708 !-----------------------------------------------------------------------------
3709       subroutine cartprint
3710
3711       use geometry_data, only: c
3712       use energy_data, only: itype
3713 !      implicit real*8 (a-h,o-z)
3714 !      include 'DIMENSIONS'
3715 !      include 'COMMON.CHAIN'
3716 !      include 'COMMON.INTERACT'
3717 !      include 'COMMON.NAMES'
3718 !      include 'COMMON.IOUNITS'
3719       integer :: i
3720
3721       write (iout,100)
3722       do i=1,nres
3723         write (iout,110) restyp(itype(i,1),1),i,c(1,i),c(2,i),&
3724           c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
3725       enddo
3726   100 format (//'              alpha-carbon coordinates       ',&
3727                 '     centroid coordinates'/ &
3728                 '       ', 6X,'X',11X,'Y',11X,'Z',&
3729                                 10X,'X',11X,'Y',11X,'Z')
3730   110 format (a,'(',i3,')',6f12.5)
3731       return
3732       end subroutine cartprint
3733 !-----------------------------------------------------------------------------
3734 !-----------------------------------------------------------------------------
3735       subroutine alloc_geo_arrays
3736 !EL Allocation of tables used by module energy
3737
3738       integer :: i,j,nres2
3739       nres2=2*nres
3740 ! commom.bounds
3741 !      common /bounds/
3742       allocate(phibound(2,nres+2)) !(2,maxres)
3743 !----------------------
3744 ! commom.chain
3745 !      common /chain/ in molread
3746 !      real(kind=8),dimension(:,:),allocatable :: c !(3,maxres2+2)
3747 !      real(kind=8),dimension(:,:),allocatable :: dc
3748       allocate(dc_old(3,0:nres2))
3749 !      if(.not.allocated(dc_norm2)) allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2)
3750       if(.not.allocated(dc_norm2)) then
3751         allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2)
3752         dc_norm2(:,:)=0.d0
3753       endif
3754 !
3755 !el      if(.not.allocated(dc_norm))
3756 !elwrite(iout,*) "jestem w alloc geo 1"
3757       if(.not.allocated(dc_norm)) then
3758         allocate(dc_norm(3,0:nres2+2)) !(3,0:maxres2)
3759         dc_norm(:,:)=0.d0
3760       endif
3761 !elwrite(iout,*) "jestem w alloc geo 1"
3762       allocate(xloc(3,nres),xrot(3,nres))
3763 !elwrite(iout,*) "jestem w alloc geo 1"
3764       xloc(:,:)=0.0D0
3765 !elwrite(iout,*) "jestem w alloc geo 1"
3766       allocate(dc_work(6*nres)) !(MAXRES6) maxres6=6*maxres
3767 !      common /rotmat/
3768       allocate(t(3,3,nres),r(3,3,nres))
3769       allocate(prod(3,3,nres),rt(3,3,nres)) !(3,3,maxres)
3770 !      common /refstruct/
3771       if(.not.allocated(cref)) allocate(cref(3,nres2+2,maxperm)) !(3,maxres2+2,maxperm)
3772 !elwrite(iout,*) "jestem w alloc geo 2"
3773       allocate(crefjlee(3,nres2+2)) !(3,maxres2+2)
3774       if(.not.allocated(chain_rep)) allocate(chain_rep(3,nres2+2,maxsym)) !(3,maxres2+2,maxsym)
3775       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3776 !      common /from_zscore/ in module.compare
3777 !----------------------
3778 ! common.local
3779 ! Inverses of the actual virtual bond lengths
3780 !      common /invlen/ in io_conf: molread or readpdb
3781 !      real(kind=8),dimension(:),allocatable :: vbld_inv !(maxres2)
3782 !----------------------
3783 ! common.var
3784 ! Store the geometric variables in the following COMMON block.
3785 !      common /var/ in readpdb or ...
3786       if(.not.allocated(theta)) allocate(theta(nres+2))
3787       if(.not.allocated(phi)) allocate(phi(nres+2))
3788       if(.not.allocated(alph)) allocate(alph(nres+2))
3789       if(.not.allocated(omeg)) allocate(omeg(nres+2))
3790       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3791       if(.not.allocated(phiref)) allocate(phiref(nres+2))
3792       if(.not.allocated(costtab)) allocate(costtab(nres))
3793       if(.not.allocated(sinttab)) allocate(sinttab(nres))
3794       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3795       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3796 !      real(kind=8),dimension(:),allocatable :: vbld !(2*maxres) in io_conf: molread or readpdb
3797       allocate(omicron(2,nres+2)) !(2,maxres)
3798       allocate(tauangle(3,nres+2)) !(3,maxres)
3799 !elwrite(iout,*) "jestem w alloc geo 3"
3800       if(.not.allocated(xxtab)) allocate(xxtab(nres))
3801       if(.not.allocated(yytab)) allocate(yytab(nres))
3802       if(.not.allocated(zztab)) allocate(zztab(nres)) !(maxres)
3803       if(.not.allocated(xxref)) allocate(xxref(nres))
3804       if(.not.allocated(yyref)) allocate(yyref(nres))
3805       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3806       allocate(ialph(nres,2)) !(maxres,2)
3807       ialph(:,1)=0
3808       ialph(:,2)=0
3809       allocate(ivar(4*nres2)) !(4*maxres2)
3810
3811 #if defined(WHAM_RUN) || defined(CLUSTER)
3812       allocate(vbld(2*nres))
3813       vbld(:)=0.d0
3814       allocate(vbld_inv(2*nres))
3815       vbld_inv(:)=0.d0
3816 #endif
3817
3818       return
3819       end subroutine alloc_geo_arrays
3820 !-----------------------------------------------------------------------------
3821 !-----------------------------------------------------------------------------
3822       subroutine returnbox
3823       integer :: allareout,i,j,k,nojumpval,chain_beg,mnum
3824       integer :: chain_end,ireturnval
3825       real*8 :: difference
3826 !C change suggested by Ana - end
3827         j=1
3828         chain_beg=1
3829 !C        do i=1,nres
3830 !C       write(*,*) 'initial', i,j,c(j,i)
3831 !C        enddo
3832 !C change suggested by Ana - begin
3833         allareout=1
3834 !C change suggested by Ana -end
3835         do i=1,nres-1
3836            mnum=molnum(i)
3837          if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
3838             .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
3839           chain_end=i
3840           if (allareout.eq.1) then
3841             ireturnval=int(c(j,i)/boxxsize)
3842             if (c(j,i).le.0) ireturnval=ireturnval-1
3843             do k=chain_beg,chain_end
3844               c(j,k)=c(j,k)-ireturnval*boxxsize
3845               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
3846             enddo
3847 !C Suggested by Ana
3848             if (chain_beg.eq.1) &
3849             dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
3850 !C Suggested by Ana -end
3851            endif
3852            chain_beg=i+1
3853            allareout=1
3854          else
3855           if (int(c(j,i)/boxxsize).eq.0) allareout=0
3856          endif
3857         enddo
3858          if (allareout.eq.1) then
3859             ireturnval=int(c(j,i)/boxxsize)
3860             if (c(j,i).le.0) ireturnval=ireturnval-1
3861             do k=chain_beg,nres
3862               c(j,k)=c(j,k)-ireturnval*boxxsize
3863               c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
3864             enddo
3865           endif
3866 !C NO JUMP
3867 !C        do i=1,nres
3868 !C        write(*,*) 'befor no jump', i,j,c(j,i)
3869 !C        enddo
3870         nojumpval=0
3871         do i=2,nres
3872            mnum=molnum(i)
3873            if (itype(i,mnum).eq.ntyp1_molec(mnum)&
3874               .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
3875              difference=abs(c(j,i-1)-c(j,i))
3876 !C             print *,'diff', difference
3877              if (difference.gt.boxxsize/2.0) then
3878                 if (c(j,i-1).gt.c(j,i)) then
3879                   nojumpval=1
3880                  else
3881                    nojumpval=-1
3882                  endif
3883               else
3884               nojumpval=0
3885               endif
3886               endif
3887               c(j,i)=c(j,i)+nojumpval*boxxsize
3888               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
3889          enddo
3890        nojumpval=0
3891         do i=2,nres
3892            mnum=molnum(i)
3893            if (itype(i,mnum).eq.ntyp1_molec(mnum) .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
3894              difference=abs(c(j,i-1)-c(j,i))
3895              if (difference.gt.boxxsize/2.0) then
3896                 if (c(j,i-1).gt.c(j,i)) then
3897                   nojumpval=1
3898                  else
3899                    nojumpval=-1
3900                  endif
3901               else
3902               nojumpval=0
3903               endif
3904              endif
3905               c(j,i)=c(j,i)+nojumpval*boxxsize
3906               c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
3907          enddo
3908
3909 !C        do i=1,nres
3910 !C        write(*,*) 'after no jump', i,j,c(j,i)
3911 !C        enddo
3912
3913 !C NOW Y dimension
3914 !C suggesed by Ana begins
3915         allareout=1
3916         j=2
3917         chain_beg=1
3918         do i=1,nres-1
3919            mnum=molnum(i)
3920          if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
3921            .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
3922           chain_end=i
3923           if (allareout.eq.1) then
3924             ireturnval=int(c(j,i)/boxysize)
3925             if (c(j,i).le.0) ireturnval=ireturnval-1
3926             do k=chain_beg,chain_end
3927               c(j,k)=c(j,k)-ireturnval*boxysize
3928              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
3929             enddo
3930 !C Suggested by Ana
3931             if (chain_beg.eq.1) &
3932             dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
3933 !C Suggested by Ana -end
3934            endif
3935            chain_beg=i+1
3936            allareout=1
3937          else
3938           if (int(c(j,i)/boxysize).eq.0) allareout=0
3939          endif
3940         enddo
3941          if (allareout.eq.1) then
3942             ireturnval=int(c(j,i)/boxysize)
3943             if (c(j,i).le.0) ireturnval=ireturnval-1
3944             do k=chain_beg,nres
3945               c(j,k)=c(j,k)-ireturnval*boxysize
3946               c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
3947             enddo
3948           endif
3949         nojumpval=0
3950         do i=2,nres
3951            mnum=molnum(i)
3952            if (itype(i,mnum).eq.ntyp1_molec(mnum)&
3953               .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
3954              difference=abs(c(j,i-1)-c(j,i))
3955              if (difference.gt.boxysize/2.0) then
3956                 if (c(j,i-1).gt.c(j,i)) then
3957                   nojumpval=1
3958                  else
3959                    nojumpval=-1
3960                  endif
3961              else
3962               nojumpval=0
3963               endif
3964            endif
3965               c(j,i)=c(j,i)+nojumpval*boxysize
3966               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
3967          enddo
3968       nojumpval=0
3969         do i=2,nres
3970            mnum=molnum(i)
3971            if (itype(i,mnum).eq.ntyp1_molec(mnum)&
3972              .and. itype(i-1,mnum).eq.ntyp1) then
3973              difference=abs(c(j,i-1)-c(j,i))
3974              if (difference.gt.boxysize/2.0) then
3975                 if (c(j,i-1).gt.c(j,i)) then
3976                   nojumpval=1
3977                  else
3978                    nojumpval=-1
3979                  endif
3980               else
3981               nojumpval=0
3982               endif
3983             endif
3984               c(j,i)=c(j,i)+nojumpval*boxysize
3985               c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
3986          enddo
3987 !C Now Z dimension
3988 !C Suggested by Ana -begins
3989         allareout=1
3990 !C Suggested by Ana -ends
3991        j=3
3992         chain_beg=1
3993         do i=1,nres-1
3994            mnum=molnum(i)
3995          if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
3996            .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
3997           chain_end=i
3998           if (allareout.eq.1) then
3999             ireturnval=int(c(j,i)/boxysize)
4000             if (c(j,i).le.0) ireturnval=ireturnval-1
4001             do k=chain_beg,chain_end
4002               c(j,k)=c(j,k)-ireturnval*boxzsize
4003               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
4004             enddo
4005 !C Suggested by Ana
4006             if (chain_beg.eq.1) dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
4007 !C Suggested by Ana -end
4008            endif
4009            chain_beg=i+1
4010            allareout=1
4011          else
4012           if (int(c(j,i)/boxzsize).eq.0) allareout=0
4013          endif
4014         enddo
4015          if (allareout.eq.1) then
4016             ireturnval=int(c(j,i)/boxzsize)
4017             if (c(j,i).le.0) ireturnval=ireturnval-1
4018             do k=chain_beg,nres
4019               c(j,k)=c(j,k)-ireturnval*boxzsize
4020               c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
4021             enddo
4022           endif
4023         nojumpval=0
4024         do i=2,nres
4025            mnum=molnum(i)
4026            if (itype(i,mnum).eq.ntyp1_molec(mnum) .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
4027              difference=abs(c(j,i-1)-c(j,i))
4028              if (difference.gt.(boxzsize/2.0)) then
4029                 if (c(j,i-1).gt.c(j,i)) then
4030                   nojumpval=1
4031                  else
4032                    nojumpval=-1
4033                  endif
4034               else
4035               nojumpval=0
4036               endif
4037             endif
4038               c(j,i)=c(j,i)+nojumpval*boxzsize
4039               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
4040          enddo
4041        nojumpval=0
4042         do i=2,nres
4043            mnum=molnum(i)
4044            if (itype(i,mnum).eq.ntyp1_molec(mnum) &
4045             .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
4046              difference=abs(c(j,i-1)-c(j,i))
4047              if (difference.gt.boxzsize/2.0) then
4048                 if (c(j,i-1).gt.c(j,i)) then
4049                   nojumpval=1
4050                  else
4051                    nojumpval=-1
4052                  endif
4053               else
4054               nojumpval=0
4055               endif
4056             endif
4057              c(j,i)=c(j,i)+nojumpval*boxzsize
4058               c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
4059          enddo
4060         do i=1,nres
4061          if (molnum(i).eq.5) then
4062           c(1,i)=dmod(c(1,i),boxxsize)
4063           if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
4064           c(2,i)=dmod(c(2,i),boxysize)
4065           if (c(2,i).lt.0) c(2,i)=c(2,i)+boxysize
4066           c(3,i)=dmod(c(3,i),boxzsize)
4067           if (c(3,i).lt.0) c(3,i)=c(3,i)+boxzsize
4068           c(1,i+nres)=dmod(c(1,i+nres),boxxsize)
4069           c(2,i+nres)=dmod(c(2,i+nres),boxysize)
4070           c(3,i+nres)=dmod(c(3,i+nres),boxzsize)
4071          endif
4072         enddo
4073         return
4074         end       subroutine returnbox
4075 !-------------------------------------------------------------------------------------------------------
4076       end module geometry