source/unres/geometry.f90

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