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