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