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