source/unres/compare.F90

   1       module compare
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use names
   5       use geometry_data
   6       use energy_data
   7       use control_data
   8 #if .not. defined WHAM_RUN && .not. defined CLUSTER
   9       use compare_data
  10       use io_base
  11       use io_config
  12       use geometry
  13       use energy
  14       use control, only: hpb_partition
  15       use minim_data
  16       use minimm, only: sc_move, minimize
  17 #endif
  18       implicit none
  19 !-----------------------------------------------------------------------------
  20 !
  21 !
  22 !-----------------------------------------------------------------------------
  23       contains
  24 #if .not. defined WHAM_RUN && .not. defined CLUSTER
  25 !-----------------------------------------------------------------------------
  26 ! contact.f
  27 !-----------------------------------------------------------------------------
  28       subroutine contact(lprint,ncont,icont,co)
  29
  30       use geometry, only:dist
  31 !      implicit real*8 (a-h,o-z)
  32 !      include 'DIMENSIONS'
  33 !      include 'COMMON.IOUNITS'
  34 !      include 'COMMON.CHAIN'
  35 !      include 'COMMON.INTERACT'
  36 !      include 'COMMON.FFIELD'
  37 !      include 'COMMON.NAMES'
  38       real(kind=8) :: facont=1.569D0  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
  39       integer :: ncont
  40       integer,dimension(2,12*nres) :: icont!(2,12*nres) !(2,maxcont)    (maxcont=12*maxres)
  41       logical :: lprint
  42 !el local variables
  43       real(kind=8) :: co,rcomp
  44       integer :: kkk,i,j,i1,i2,it1,it2,iti,itj
  45
  46       ncont=0
  47       kkk=3
  48       do i=nnt+kkk,nct
  49         iti=iabs(itype(i,1))
  50         do j=nnt,i-kkk
  51           itj=iabs(itype(j,1))
  52           if (ipot.ne.4) then
  53 !           rcomp=sigmaii(iti,itj)+1.0D0
  54             rcomp=facont*sigmaii(iti,itj)
  55           else
  56 !           rcomp=sigma(iti,itj)+1.0D0
  57             rcomp=facont*sigma(iti,itj)
  58           endif
  59 !         rcomp=6.5D0
  60 !         print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
  61           if (dist(nres+i,nres+j).lt.rcomp) then
  62             ncont=ncont+1
  63             icont(1,ncont)=i
  64             icont(2,ncont)=j
  65           endif
  66         enddo
  67       enddo
  68       if (lprint) then
  69         write (iout,'(a)') 'Contact map:'
  70         do i=1,ncont
  71           i1=icont(1,i)
  72           i2=icont(2,i)
  73           it1=itype(i1,1)
  74           it2=itype(i2,1)
  75           write (iout,'(i3,2x,a,i4,2x,a,i4)') &
  76            i,restyp(it1,1),i1,restyp(it2,1),i2
  77         enddo
  78       endif
  79       co = 0.0d0
  80       do i=1,ncont
  81         co = co + dfloat(iabs(icont(1,i)-icont(2,i)))
  82       enddo
  83       co = co / (nres*ncont)
  84       return
  85       end subroutine contact
  86 !-----------------------------------------------------------------------------
  87       real(kind=8) function contact_fract(ncont,ncont_ref,icont,icont_ref)
  88
  89 !      implicit real*8 (a-h,o-z)
  90 !      include 'DIMENSIONS'
  91 !      include 'COMMON.IOUNITS'
  92       integer :: ncont,ncont_ref
  93       integer,dimension(2,12*nres) :: icont,icont_ref   !(2,12*nres) (2,maxcont)        (maxcont=12*maxres)
  94 !el local variables
  95       integer :: i,j,nmatch
  96       nmatch=0
  97 !     print *,'ncont=',ncont,' ncont_ref=',ncont_ref
  98 !     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
  99 !     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
 100 !     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
 101 !     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
 102       do i=1,ncont
 103         do j=1,ncont_ref
 104           if (icont(1,i).eq.icont_ref(1,j) .and. &
 105               icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
 106         enddo
 107       enddo
 108 !     print *,' nmatch=',nmatch
 109 !     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
 110       contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
 111       return
 112       end function contact_fract
 113 !-----------------------------------------------------------------------------
 114       real(kind=8) function contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
 115
 116 !      implicit real*8 (a-h,o-z)
 117 !      include 'DIMENSIONS'
 118 !      include 'COMMON.IOUNITS'
 119       integer :: ncont,ncont_ref
 120       integer,dimension(2,12*nres) :: icont,icont_ref   !(2,12*nres) (2,maxcont)        (maxcont=12*maxres)
 121 !el local variables
 122       integer :: i,j,nmatch
 123       nmatch=0
 124 !     print *,'ncont=',ncont,' ncont_ref=',ncont_ref
 125 !     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
 126 !     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
 127 !     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
 128 !     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
 129       do i=1,ncont
 130         do j=1,ncont_ref
 131           if (icont(1,i).eq.icont_ref(1,j) .and. &
 132               icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
 133         enddo
 134       enddo
 135 !     print *,' nmatch=',nmatch
 136 !     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
 137       contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont)
 138       return
 139       end function contact_fract_nn
 140 !-----------------------------------------------------------------------------
 141       subroutine hairpin(lprint,nharp,iharp)
 142
 143       use geometry, only:dist
 144 !      implicit real*8 (a-h,o-z)
 145 !      include 'DIMENSIONS'
 146 !      include 'COMMON.IOUNITS'
 147 !      include 'COMMON.CHAIN'
 148 !      include 'COMMON.INTERACT'
 149 !      include 'COMMON.FFIELD'
 150 !      include 'COMMON.NAMES'
 151       integer :: ncont
 152       integer,dimension(2,12*nres) :: icont     !(2,maxcont)    (maxcont=12*maxres)
 153       integer :: nharp
 154       integer,dimension(4,nres/3) :: iharp      !(4,nres/3)(4,maxres/3)
 155       logical :: lprint,not_done
 156       real(kind=8) :: rcomp=6.0d0
 157 !el local variables
 158       integer :: i,j,kkk,k,i1,i2,it1,it2,j1,ii1,jj1
 159 !      allocate(icont(2,12*nres))
 160
 161       ncont=0
 162       kkk=0
 163 !     print *,'nnt=',nnt,' nct=',nct
 164       do i=nnt,nct-3
 165         do k=1,3
 166           c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
 167         enddo
 168         do j=i+2,nct-1
 169           do k=1,3
 170             c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
 171           enddo
 172           if (dist(2*nres+1,2*nres+2).lt.rcomp) then
 173             ncont=ncont+1
 174             icont(1,ncont)=i
 175             icont(2,ncont)=j
 176           endif
 177         enddo
 178       enddo
 179       if (lprint) then
 180         write (iout,'(a)') 'PP contact map:'
 181         do i=1,ncont
 182           i1=icont(1,i)
 183           i2=icont(2,i)
 184           it1=itype(i1,1)
 185           it2=itype(i2,1)
 186           write (iout,'(i3,2x,a,i4,2x,a,i4)') &
 187            i,restyp(it1,1),i1,restyp(it2,1),i2
 188         enddo
 189       endif
 190 ! finding hairpins
 191       nharp=0
 192       do i=1,ncont
 193         i1=icont(1,i)
 194         j1=icont(2,i)
 195         if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then
 196 !          write (iout,*) "found turn at ",i1,j1
 197           ii1=i1
 198           jj1=j1
 199           not_done=.true.
 200           do while (not_done)
 201             i1=i1-1
 202             j1=j1+1
 203             do j=1,ncont
 204               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
 205             enddo
 206             not_done=.false.
 207   10        continue
 208 !            write (iout,*) i1,j1,not_done
 209           enddo
 210           i1=i1+1
 211           j1=j1-1
 212           if (j1-i1.gt.4) then
 213             nharp=nharp+1
 214             iharp(1,nharp)=i1
 215             iharp(2,nharp)=j1
 216             iharp(3,nharp)=ii1
 217             iharp(4,nharp)=jj1
 218 !            write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4)
 219           endif
 220         endif
 221       enddo
 222 !      do i=1,nharp
 223 !            write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
 224 !      enddo
 225       if (lprint) then
 226       write (iout,*) "Hairpins:"
 227       do i=1,nharp
 228         i1=iharp(1,i)
 229         j1=iharp(2,i)
 230         ii1=iharp(3,i)
 231         jj1=iharp(4,i)
 232         write (iout,*)
 233         write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=i1,ii1)
 234         write (iout,'(20(a,i3,1x))') (restyp(itype(k,1),1),k,k=j1,jj1,-1)
 235 !        do k=jj1,j1,-1
 236 !         write (iout,'(a,i3,$)') restyp(itype(k,1)),k
 237 !        enddo
 238       enddo
 239       endif
 240       return
 241       end subroutine hairpin
 242 !-----------------------------------------------------------------------------
 243 ! elecont.f
 244 !-----------------------------------------------------------------------------
 245       subroutine elecont(lprint,ncont,icont)
 246
 247 !      implicit real*8 (a-h,o-z)
 248 !      include 'DIMENSIONS'
 249 !      include 'COMMON.IOUNITS'
 250 !      include 'COMMON.CHAIN'
 251 !      include 'COMMON.INTERACT'
 252 !      include 'COMMON.LOCAL'
 253 !      include 'COMMON.FFIELD'
 254 !      include 'COMMON.NAMES'
 255       logical :: lprint
 256       real(kind=8),dimension(2,2) :: elpp_6,elpp_3,ael6_,ael3_
 257       real(kind=8) :: ael6_i,ael3_i
 258       real(kind=8),dimension(2,2) :: app_,bpp_,rpp_
 259       integer :: ncont
 260       integer,dimension(2,12*nres) :: icont     !(2,12*nres)(2,maxcont) (maxcont=12*maxres)
 261       real(kind=8),dimension(12*nres) :: econt  !(maxcont)
 262 !el local variables
 263       integer :: i,j,k,iteli,itelj,i1,i2,it1,it2,ic1,ic2
 264       real(kind=8) :: elcutoff,elecutoff_14,rri,ees,evdw
 265       real(kind=8) :: xi,yi,zi,dxi,dyi,dzi,aaa,bbb
 266       real(kind=8) :: xmedi,ymedi,zmedi
 267       real(kind=8) :: xj,yj,zj,dxj,dyj,dzj,rrmij,rmij,r3ij,r6ij
 268       real(kind=8) :: vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,&
 269                  evdwij,el1,el2,eesij,ene
 270 !
 271 ! Load the constants of peptide bond - peptide bond interactions.
 272 ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
 273 ! proline) - determined by averaging ECEPP energy.
 274 !
 275 ! as of 7/06/91.
 276 !
 277 !      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
 278       data rpp_    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
 279       data elpp_6  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
 280       data elpp_3  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
 281
 282 !el      allocate(econt(12*nres))       !(maxcont)
 283
 284       elcutoff = -0.3d0
 285       elecutoff_14 = -0.5d0
 286       if (lprint) write (iout,'(a)') &
 287         "Constants of electrostatic interaction energy expression."
 288       do i=1,2
 289         do j=1,2
 290         rri=rpp_(i,j)**6
 291         app_(i,j)=epp(i,j)*rri*rri
 292         bpp_(i,j)=-2.0*epp(i,j)*rri
 293         ael6_(i,j)=elpp_6(i,j)*4.2**6
 294         ael3_(i,j)=elpp_3(i,j)*4.2**3
 295         if (lprint) &
 296         write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),&
 297                                      ael3_(i,j)
 298         enddo
 299       enddo
 300       ncont=0
 301       ees=0.0
 302       evdw=0.0
 303       print *, "nntt,nct",nnt,nct-2
 304       do 1 i=nnt,nct-2
 305         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
 306         xi=c(1,i)
 307         yi=c(2,i)
 308         zi=c(3,i)
 309         dxi=c(1,i+1)-c(1,i)
 310         dyi=c(2,i+1)-c(2,i)
 311         dzi=c(3,i+1)-c(3,i)
 312         xmedi=xi+0.5*dxi
 313         ymedi=yi+0.5*dyi
 314         zmedi=zi+0.5*dzi
 315         do 4 j=i+2,nct-1
 316           if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) goto 4
 317           iteli=itel(i)
 318           itelj=itel(j)
 319           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
 320           if (iteli.eq.2 .and. itelj.eq.2) goto 4
 321           aaa=app_(iteli,itelj)
 322           bbb=bpp_(iteli,itelj)
 323           ael6_i=ael6_(iteli,itelj)
 324           ael3_i=ael3_(iteli,itelj)
 325           dxj=c(1,j+1)-c(1,j)
 326           dyj=c(2,j+1)-c(2,j)
 327           dzj=c(3,j+1)-c(3,j)
 328           xj=c(1,j)+0.5*dxj-xmedi
 329           yj=c(2,j)+0.5*dyj-ymedi
 330           zj=c(3,j)+0.5*dzj-zmedi
 331           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
 332           rmij=sqrt(rrmij)
 333           r3ij=rrmij*rmij
 334           r6ij=r3ij*r3ij
 335           vrmij=vblinv*rmij
 336           cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
 337           cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
 338           cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
 339           fac=cosa-3.0*cosb*cosg
 340           ev1=aaa*r6ij*r6ij
 341           ev2=bbb*r6ij
 342           fac3=ael6_i*r6ij
 343           fac4=ael3_i*r3ij
 344           evdwij=ev1+ev2
 345           el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
 346           el2=fac4*fac
 347           eesij=el1+el2
 348           if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
 349               j.eq.i+2 .and. eesij.le.elecutoff_14) then
 350              ncont=ncont+1
 351              icont(1,ncont)=i
 352              icont(2,ncont)=j
 353              econt(ncont)=eesij
 354           endif
 355           ees=ees+eesij
 356           evdw=evdw+evdwij
 357     4   continue
 358     1 continue
 359       if (lprint) then
 360         write (iout,*) 'Total average electrostatic energy: ',ees
 361         write (iout,*) 'VDW energy between peptide-group centers: ',evdw
 362         write (iout,*)
 363         write (iout,*) 'Electrostatic contacts before pruning: '
 364         do i=1,ncont
 365           i1=icont(1,i)
 366           i2=icont(2,i)
 367           it1=itype(i1,1)
 368           it2=itype(i2,1)
 369           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
 370            i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
 371         enddo
 372       endif
 373 ! For given residues keep only the contacts with the greatest energy.
 374       i=0
 375       do while (i.lt.ncont)
 376         i=i+1
 377         ene=econt(i)
 378         ic1=icont(1,i)
 379         ic2=icont(2,i)
 380         j=i
 381         do while (j.lt.ncont)
 382           j=j+1
 383           if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
 384               ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
 385 !            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
 386 !     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
 387             if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
 388               if (ic1.eq.icont(1,j)) then
 389                 do k=1,ncont
 390                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j) &
 391                      .and. iabs(icont(1,k)-ic1).le.2 .and. &
 392                      econt(k).lt.econt(j) ) goto 21
 393                 enddo
 394               else if (ic2.eq.icont(2,j) ) then
 395                 do k=1,ncont
 396                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j) &
 397                      .and. iabs(icont(2,k)-ic2).le.2 .and. &
 398                      econt(k).lt.econt(j) ) goto 21
 399                 enddo
 400               endif
 401 ! Remove ith contact
 402               do k=i+1,ncont
 403                 icont(1,k-1)=icont(1,k)
 404                 icont(2,k-1)=icont(2,k)
 405                 econt(k-1)=econt(k)
 406               enddo
 407               i=i-1
 408               ncont=ncont-1
 409 !              write (iout,*) "ncont",ncont
 410 !              do k=1,ncont
 411 !                write (iout,*) icont(1,k),icont(2,k)
 412 !              enddo
 413               goto 20
 414             else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
 415             then
 416               if (ic1.eq.icont(1,j)) then
 417                 do k=1,ncont
 418                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
 419                      .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
 420                      econt(k).lt.econt(i) ) goto 21
 421                 enddo
 422               else if (ic2.eq.icont(2,j) ) then
 423                 do k=1,ncont
 424                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
 425                      .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
 426                      econt(k).lt.econt(i) ) goto 21
 427                 enddo
 428               endif
 429 ! Remove jth contact
 430               do k=j+1,ncont
 431                 icont(1,k-1)=icont(1,k)
 432                 icont(2,k-1)=icont(2,k)
 433                 econt(k-1)=econt(k)
 434               enddo
 435               ncont=ncont-1
 436 !              write (iout,*) "ncont",ncont
 437 !              do k=1,ncont
 438 !                write (iout,*) icont(1,k),icont(2,k)
 439 !              enddo
 440               j=j-1
 441             endif
 442           endif
 443    21     continue
 444         enddo
 445    20   continue
 446       enddo
 447       if (lprint) then
 448         write (iout,*)
 449         write (iout,*) 'Electrostatic contacts after pruning: '
 450         do i=1,ncont
 451           i1=icont(1,i)
 452           i2=icont(2,i)
 453           it1=itype(i1,1)
 454           it2=itype(i2,1)
 455           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
 456            i,restyp(it1,1),i1,restyp(it2,1),i2,econt(i)
 457         enddo
 458       endif
 459       return
 460       end subroutine elecont
 461 !-----------------------------------------------------------------------------
 462       subroutine secondary2(lprint)
 463
 464 !      implicit real*8 (a-h,o-z)
 465 !      include 'DIMENSIONS'
 466 !      include 'COMMON.CHAIN'
 467 !      include 'COMMON.IOUNITS'
 468 !      include 'COMMON.DISTFIT'
 469 !      include 'COMMON.VAR'
 470 !      include 'COMMON.GEO'
 471 !      include 'COMMON.CONTROL'
 472       integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,&
 473              iii1,jjj1
 474       integer,dimension(2,12*nres) :: icont     !(2,maxcont)    (maxcont=12*maxres)
 475       integer,dimension(nres,0:4) :: isec       !(maxres,4)
 476       integer,dimension(nres) :: nsec   !(maxres)
 477       logical :: lprint,not_done        !,freeres
 478       real(kind=8) :: p1,p2
 479 !el      external freeres
 480
 481 !el      allocate(icont(2,12*nres),isec(nres,4),nsec(nres))
 482
 483       if(.not.dccart) call chainbuild_cart
 484       if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
 485 !d      call write_pdb(99,'sec structure',0d0)
 486       ncont=0
 487       nbfrag=0
 488       nhfrag=0
 489       do i=1,nres
 490         isec(i,1)=0
 491         isec(i,2)=0
 492         nsec(i)=0
 493       enddo
 494
 495       call elecont(lprint,ncont,icont)
 496       print *,"after elecont"
 497       if (nres_molec(1).eq.0) return
 498
 499 ! finding parallel beta
 500 !d      write (iout,*) '------- looking for parallel beta -----------'
 501       nbeta=0
 502       nstrand=0
 503       do i=1,ncont
 504         i1=icont(1,i)
 505         j1=icont(2,i)
 506         if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
 507           ii1=i1
 508           jj1=j1
 509 !d          write (iout,*) i1,j1
 510           not_done=.true.
 511           do while (not_done)
 512            i1=i1+1
 513            j1=j1+1
 514             do j=1,ncont
 515               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
 516                    freeres(i1,j1,nsec,isec)) goto 5
 517             enddo
 518             not_done=.false.
 519   5         continue
 520 !d            write (iout,*) i1,j1,not_done
 521           enddo
 522           j1=j1-1
 523           i1=i1-1
 524           if (i1-ii1.gt.1) then
 525             ii1=max0(ii1-1,1)
 526             jj1=max0(jj1-1,1)
 527             nbeta=nbeta+1
 528             if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
 529                      nbeta,ii1,i1,jj1,j1
 530
 531             nbfrag=nbfrag+1
 532             bfrag(1,nbfrag)=ii1+1
 533             bfrag(2,nbfrag)=i1+1
 534             bfrag(3,nbfrag)=jj1+1
 535             bfrag(4,nbfrag)=min0(j1+1,nres)
 536
 537             do ij=ii1,i1
 538              nsec(ij)=nsec(ij)+1
 539              isec(ij,nsec(ij))=nbeta
 540             enddo
 541             do ij=jj1,j1
 542              nsec(ij)=nsec(ij)+1
 543              isec(ij,nsec(ij))=nbeta
 544             enddo
 545
 546            if(lprint) then
 547             nstrand=nstrand+1
 548             if (nbeta.le.9) then
 549               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
 550                 "DefPropRes 'strand",nstrand,&
 551                 "' 'num = ",ii1-1,"..",i1-1,"'"
 552             else
 553               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
 554                 "DefPropRes 'strand",nstrand,&
 555                 "' 'num = ",ii1-1,"..",i1-1,"'"
 556             endif
 557             nstrand=nstrand+1
 558             if (nbeta.le.9) then
 559               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
 560                 "DefPropRes 'strand",nstrand,&
 561                 "' 'num = ",jj1-1,"..",j1-1,"'"
 562             else
 563               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
 564                 "DefPropRes 'strand",nstrand,&
 565                 "' 'num = ",jj1-1,"..",j1-1,"'"
 566             endif
 567               write(12,'(a8,4i4)') &
 568                 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
 569            endif
 570           endif
 571         endif
 572       enddo
 573
 574 ! finding alpha or 310 helix
 575       nhelix=0
 576       do i=1,ncont
 577         i1=icont(1,i)
 578         j1=icont(2,i)
 579         p1=phi(i1+2)*rad2deg
 580         p2=0.0
 581         if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
 582
 583
 584         if (j1.eq.i1+3 .and. &
 585              ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
 586              ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
 587 !d          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
 588 !o          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
 589           ii1=i1
 590           jj1=j1
 591           if (nsec(ii1).eq.0) then
 592             not_done=.true.
 593           else
 594             not_done=.false.
 595           endif
 596           do while (not_done)
 597             i1=i1+1
 598             j1=j1+1
 599             do j=1,ncont
 600               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
 601             enddo
 602             not_done=.false.
 603   10        continue
 604             p1=phi(i1+2)*rad2deg
 605             p2=phi(j1+2)*rad2deg
 606             if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
 607                                     not_done=.false.
 608 !d
 609           enddo
 610           j1=j1+1
 611           if (j1-ii1.gt.5) then
 612             nhelix=nhelix+1
 613 !d
 614
 615             nhfrag=nhfrag+1
 616             hfrag(1,nhfrag)=ii1
 617             hfrag(2,nhfrag)=j1
 618
 619             do ij=ii1,j1
 620              nsec(ij)=-1
 621             enddo
 622            if (lprint) then
 623             write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
 624             if (nhelix.le.9) then
 625               write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
 626                 "DefPropRes 'helix",nhelix,&
 627                 "' 'num = ",ii1-1,"..",j1-2,"'"
 628             else
 629               write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
 630                 "DefPropRes 'helix",nhelix,&
 631                 "' 'num = ",ii1-1,"..",j1-2,"'"
 632             endif
 633            endif
 634           endif
 635         endif
 636       enddo
 637       if (nhelix.gt.0.and.lprint) then
 638         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
 639         do i=2,nhelix
 640          if (nhelix.le.9) then
 641           write(12,'(a8,i1,$)') " | helix",i
 642          else
 643           write(12,'(a8,i2,$)') " | helix",i
 644          endif
 645         enddo
 646         write(12,'(a1)') "'"
 647       endif
 648
 649
 650 ! finding antiparallel beta
 651 !d      write (iout,*) '--------- looking for antiparallel beta ---------'
 652
 653       do i=1,ncont
 654         i1=icont(1,i)
 655         j1=icont(2,i)
 656         if (freeres(i1,j1,nsec,isec)) then
 657           ii1=i1
 658           jj1=j1
 659 !d          write (iout,*) i1,j1
 660
 661           not_done=.true.
 662           do while (not_done)
 663            i1=i1+1
 664            j1=j1-1
 665             do j=1,ncont
 666               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
 667                    freeres(i1,j1,nsec,isec)) goto 6
 668             enddo
 669             not_done=.false.
 670   6         continue
 671 !d            write (iout,*) i1,j1,not_done
 672           enddo
 673           i1=i1-1
 674           j1=j1+1
 675           if (i1-ii1.gt.1) then
 676
 677             nbfrag=nbfrag+1
 678             bfrag(1,nbfrag)=ii1
 679             bfrag(2,nbfrag)=min0(i1+1,nres)
 680             bfrag(3,nbfrag)=min0(jj1+1,nres)
 681             bfrag(4,nbfrag)=j1
 682
 683             nbeta=nbeta+1
 684             iii1=max0(ii1-1,1)
 685             do ij=iii1,i1
 686              nsec(ij)=nsec(ij)+1
 687              if (nsec(ij).le.2) then
 688               isec(ij,nsec(ij))=nbeta
 689              endif
 690             enddo
 691             jjj1=max0(j1-1,1)
 692             do ij=jjj1,jj1
 693              nsec(ij)=nsec(ij)+1
 694              if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
 695               isec(ij,nsec(ij))=nbeta
 696              endif
 697             enddo
 698
 699
 700            if (lprint) then
 701             write (iout,'(a,i3,4i4)')'antiparallel beta',&
 702                          nbeta,ii1-1,i1,jj1,j1-1
 703             nstrand=nstrand+1
 704             if (nstrand.le.9) then
 705               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
 706                 "DefPropRes 'strand",nstrand,&
 707                 "' 'num = ",ii1-2,"..",i1-1,"'"
 708             else
 709               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
 710                 "DefPropRes 'strand",nstrand,&
 711                 "' 'num = ",ii1-2,"..",i1-1,"'"
 712             endif
 713             nstrand=nstrand+1
 714             if (nstrand.le.9) then
 715               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
 716                 "DefPropRes 'strand",nstrand,&
 717                 "' 'num = ",j1-2,"..",jj1-1,"'"
 718             else
 719               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
 720                 "DefPropRes 'strand",nstrand,&
 721                 "' 'num = ",j1-2,"..",jj1-1,"'"
 722             endif
 723               write(12,'(a8,4i4)') &
 724                 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
 725            endif
 726           endif
 727         endif
 728       enddo
 729
 730       if (nstrand.gt.0.and.lprint) then
 731         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
 732         do i=2,nstrand
 733          if (i.le.9) then
 734           write(12,'(a9,i1,$)') " | strand",i
 735          else
 736           write(12,'(a9,i2,$)') " | strand",i
 737          endif
 738         enddo
 739         write(12,'(a1)') "'"
 740       endif
 741
 742
 743
 744       if (lprint) then
 745        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
 746        write(12,'(a20)') "XMacStand ribbon.mac"
 747
 748       if (nres_molec(1).eq.0) return
 749        write(iout,*) 'UNRES seq:'
 750        do j=1,nbfrag
 751         write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
 752        enddo
 753
 754        do j=1,nhfrag
 755         write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
 756        enddo
 757       endif
 758
 759       return
 760       end subroutine secondary2
 761 #endif
 762 !-----------------------------------------------------------------------------
 763       logical function freeres(i,j,nsec,isec)
 764
 765 !      implicit real*8 (a-h,o-z)
 766 !      include 'DIMENSIONS'
 767       integer,dimension(nres,4) :: isec         !(maxres,4)
 768       integer,dimension(nres) :: nsec           !(maxres)
 769
 770 !el local variables
 771       integer :: i,j,k,l
 772
 773       freeres=.false.
 774 #ifndef WHAM_RUN
 775       if (nsec(i).lt.0.or.nsec(j).lt.0) return
 776 #endif
 777       if (nsec(i).gt.1.or.nsec(j).gt.1) return
 778       do k=1,nsec(i)
 779         do l=1,nsec(j)
 780           if (isec(i,k).eq.isec(j,l)) return
 781         enddo
 782       enddo
 783       freeres=.true.
 784       return
 785       end function freeres
 786 !-----------------------------------------------------------------------------
 787 ! readrtns_CSA.F
 788 !-----------------------------------------------------------------------------
 789       logical function seq_comp(itypea,itypeb,length)
 790
 791 !el      implicit none
 792       integer :: length,itypea(length),itypeb(length)
 793       integer :: i
 794       do i=1,length
 795         if (itypea(i).ne.itypeb(i)) then
 796           seq_comp=.false.
 797           return
 798         endif
 799       enddo
 800       seq_comp=.true.
 801       return
 802       end function seq_comp
 803 #ifndef WHAM_RUN
 804 !-----------------------------------------------------------------------------
 805 ! rmsd.F
 806 !-----------------------------------------------------------------------------
 807       subroutine rms_nac_nnc(rms,frac,frac_nn,co,lprn)
 808
 809 !        implicit real*8 (a-h,o-z)
 810 !        include 'DIMENSIONS'
 811 !        include 'COMMON.CHAIN'
 812 !        include 'COMMON.CONTACTS'
 813 !        include 'COMMON.IOUNITS'
 814         real(kind=8) :: przes(3),obr(3,3)
 815         logical :: non_conv,lprn
 816         real(kind=8) :: rms,frac,frac_nn,co
 817 !        call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
 818 !     &             obr,non_conv)
 819 !        rms=dsqrt(rms)
 820         call rmsd(rms)
 821 !        print *,"before contact"
 822 !elte(iout,*) "rms_nacc before contact"
 823         call contact(.false.,ncont,icont,co)
 824         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 825         frac_nn=contact_fract_nn(ncont,ncont_ref,icont,icont_ref)
 826         if (lprn) write (iout,'(a,f8.3/a,f8.3/a,f8.3/a,f8.3)') &
 827           'RMS deviation from the reference structure:',rms,&
 828           ' % of native contacts:',frac*100,&
 829           ' % of nonnative contacts:',frac_nn*100,&
 830           ' contact order:',co
 831
 832       return
 833       end subroutine rms_nac_nnc
 834 !-----------------------------------------------------------------------------
 835       subroutine rmsd(drms)
 836
 837       use regularize_, only:fitsq
 838 !      implicit real*8 (a-h,o-z)
 839 !      include 'DIMENSIONS'
 840 #ifdef MPI
 841       include 'mpif.h'
 842 #endif
 843 !      include 'COMMON.CHAIN'
 844 !      include 'COMMON.IOUNITS'
 845 !      include 'COMMON.INTERACT'
 846 !      include 'COMMON.CONTROL'
 847       logical :: non_conv
 848       real(kind=8) :: przes(3),obrot(3,3)
 849       real(kind=8),dimension(3,2*nres+2) :: ccopy,crefcopy      !(3,maxres2+2) maxres2=2*maxres
 850
 851 !el local variables
 852       real(kind=8) :: drms,rminroz,roznica
 853       integer :: i,j,iatom,kkk,iti,k
 854
 855 !el      allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2))       !(3,maxres2+2) maxres2=2*maxres
 856
 857       nperm=1
 858       do i=1,symetr
 859       nperm=nperm*i
 860       enddo
 861       iatom=0
 862       rminroz=100d2
 863 !      print *,"nz_start",nz_start," nz_end",nz_end
 864 !      if (symetr.le.1) then
 865       do kkk=1,nperm
 866 !      do i=nz_start,nz_end
 867 !        iatom=iatom+1
 868 !        iti=itype(i,1)
 869 !        do k=1,3
 870 !         ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
 871 !         crefcopy(k,iatom,kkk)=cref(k,i,kkk)
 872 !        enddo
 873 !        if (iz_sc.eq.1.and.iti.ne.10) then
 874 !          iatom=iatom+1
 875 !          do k=1,3
 876 !           ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
 877 !           crefcopy(k,iatom,kkk)=cref(k,nres+i,kkk)
 878 !          enddo
 879 !        endif
 880 !      enddo
 881 !      else
 882 !      do kkk=1,nperm
 883       iatom=0
 884       print *,nz_start,nz_end,nstart_seq-nstart_sup
 885       do i=nz_start,nz_end
 886         iatom=iatom+1
 887         iti=itype(i,1)
 888         do k=1,3
 889          ccopy(k,iatom)=c(k,i+nstart_seq-nstart_sup)
 890          crefcopy(k,iatom)=cref(k,i,kkk)
 891         enddo
 892         if (iz_sc.eq.1.and.iti.ne.10) then
 893           iatom=iatom+1
 894           do k=1,3
 895            ccopy(k,iatom)=c(k,nres+i+nstart_seq-nstart_sup)
 896            crefcopy(k,iatom)=cref(k,nres+i,kkk)
 897           enddo
 898         endif
 899       enddo
 900 !      enddo
 901 !      endif
 902
 903 ! ----- diagnostics
 904 !         do kkk=1,nperm
 905 !          write (iout,*) 'Ccopy and CREFcopy'
 906 !          print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),
 907 !     &           (crefcopy(j,k),j=1,3),k=1,iatom)
 908 !          write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),
 909 !     &           (crefcopy(j,k),j=1,3),k=1,iatom)
 910 !         enddo
 911 ! ----- end diagnostics
 912 !      do kkk=1,nperm
 913       call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
 914                                             przes,obrot,non_conv)
 915       if (non_conv) then
 916           print *,'Problems in FITSQ!!! rmsd'
 917           write (iout,*) 'Problems in FITSQ!!! rmsd'
 918           print *,'Ccopy and CREFcopy'
 919           write (iout,*) 'Ccopy and CREFcopy'
 920           print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
 921                  (crefcopy(j,k),j=1,3),k=1,iatom)
 922           write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
 923                  (crefcopy(j,k),j=1,3),k=1,iatom)
 924 #ifdef MPI
 925 !          call mpi_abort(mpi_comm_world,ierror,ierrcode)
 926            roznica=100.0d10
 927 #else
 928           stop
 929 #endif
 930        endif
 931 !       write (iout,*) "roznica", roznica,kkk
 932        if (roznica.le.rminroz) rminroz=roznica
 933        enddo
 934        drms=dsqrt(dabs(rminroz))
 935 ! ---- diagnostics
 936 !        write (iout,*) "nperm,symetr", nperm,symetr
 937 ! ---- end diagnostics
 938       return
 939       end subroutine rmsd
 940 !-----------------------------------------------------------------------------
 941       subroutine rmsd_csa(drms)
 942
 943       use regularize_, only:fitsq
 944 !      implicit real*8 (a-h,o-z)
 945 !      include 'DIMENSIONS'
 946 #ifdef MPI
 947       include 'mpif.h'
 948 #endif
 949 !      include 'COMMON.CHAIN'
 950 !      include 'COMMON.IOUNITS'
 951 !      include 'COMMON.INTERACT'
 952       logical :: non_conv
 953       real(kind=8) :: przes(3),obrot(3,3)
 954       real(kind=8),dimension(:,:),allocatable :: ccopy,crefcopy !(3,maxres2+2) maxres2=2*maxres
 955       integer :: kkk,iatom,ierror,ierrcode
 956
 957 !el local variables
 958       integer ::i,j,k,iti
 959       real(kind=8) :: drms,roznica
 960
 961       allocate(ccopy(3,2*nres+2),crefcopy(3,2*nres+2))  !(3,maxres2+2) maxres2=2*maxres
 962
 963       kkk=1
 964       iatom=0
 965       do i=nz_start,nz_end
 966         iatom=iatom+1
 967         iti=itype(i,1)
 968         do k=1,3
 969          ccopy(k,iatom)=c(k,i)
 970          crefcopy(k,iatom)=crefjlee(k,i)
 971         enddo
 972         if (iz_sc.eq.1.and.iti.ne.10) then
 973           iatom=iatom+1
 974           do k=1,3
 975            ccopy(k,iatom)=c(k,nres+i)
 976            crefcopy(k,iatom)=crefjlee(k,nres+i)
 977           enddo
 978         endif
 979       enddo
 980
 981       call fitsq(roznica,ccopy(1,1),crefcopy(1,1),iatom,&
 982                                             przes,obrot,non_conv)
 983       if (non_conv) then
 984           print *,'Problems in FITSQ!!! rmsd_csa'
 985           write (iout,*) 'Problems in FITSQ!!! rmsd_csa'
 986           print *,'Ccopy and CREFcopy'
 987           write (iout,*) 'Ccopy and CREFcopy'
 988           print '(i5,3f10.5,5x,3f10.5)',(k,(ccopy(j,k),j=1,3),&
 989                  (crefcopy(j,k),j=1,3),k=1,iatom)
 990           write (iout,'(i5,3f10.5,5x,3f10.5)') (k,(ccopy(j,k),j=1,3),&
 991                  (crefcopy(j,k),j=1,3),k=1,iatom)
 992 #ifdef MPI
 993           call mpi_abort(mpi_comm_world,ierror,ierrcode)
 994 #else
 995           stop
 996 #endif
 997        endif
 998        drms=dsqrt(dabs(roznica))
 999       return
1000       end subroutine rmsd_csa
1001 !-----------------------------------------------------------------------------
1002 ! test.F
1003 !-----------------------------------------------------------------------------
1004       subroutine test
1005
1006 !el      use minim
1007       use geometry, only:pinorm
1008       use random, only:ran_number,iran_num
1009 !      implicit real*8 (a-h,o-z)
1010 !      include 'DIMENSIONS'
1011       include 'mpif.h'
1012 !      include 'COMMON.GEO'
1013 !      include 'COMMON.VAR'
1014 !      include 'COMMON.INTERACT'
1015 !      include 'COMMON.IOUNITS'
1016 !      include 'COMMON.DISTFIT'
1017 !      include 'COMMON.SBRIDGE'
1018 !      include 'COMMON.CONTROL'
1019 !      include 'COMMON.FFIELD'
1020 !      include 'COMMON.MINIM'
1021 !      include 'COMMON.CHAIN'
1022       real(kind=8) :: time0,time1
1023       real(kind=8) :: energy(0:n_ene),ee
1024       real(kind=8),dimension(6*nres) :: var,var1        !(maxvar) (maxvar=6*maxres)
1025       integer :: j1,j2,jr,i,iretcode,nfun,nft_sc
1026       logical :: debug,accepted
1027       real(kind=8) :: etot,rms,da,temp,betbol,etot0,d,phiold,&
1028                  xxr,xxh
1029       debug=.true.
1030 !el      allocate(var(6*nres),var1(6*nres))     !(maxvar) (maxvar=6*maxres)
1031
1032       call geom_to_var(nvar,var1)
1033       call chainbuild
1034       call etotal(energy)
1035       etot=energy(0)
1036       call rmsd(rms)
1037       write(iout,*) 'etot=',0,etot,rms
1038       call secondary2(.false.)
1039
1040       call write_pdb(0,'first structure',etot)
1041
1042       j1=13
1043       j2=21
1044       da=180.0*deg2rad
1045
1046
1047
1048        temp=3000.0d0
1049        betbol=1.0D0/(1.9858D-3*temp)
1050        jr=iran_num(j1,j2)
1051        d=ran_number(-pi,pi)
1052 !       phi(jr)=pinorm(phi(jr)+d)
1053        call chainbuild
1054        call etotal(energy)
1055        etot0=energy(0)
1056        call rmsd(rms)
1057        write(iout,*) 'etot=',1,etot0,rms
1058        call write_pdb(1,'perturb structure',etot0)
1059
1060       do i=2,500,2
1061        jr=iran_num(j1,j2)
1062        d=ran_number(-da,da)
1063        phiold=phi(jr)
1064        phi(jr)=pinorm(phi(jr)+d)
1065        call chainbuild
1066        call etotal(energy)
1067        etot=energy(0)
1068
1069        if (etot.lt.etot0) then
1070           accepted=.true.
1071        else
1072           accepted=.false.
1073           xxr=ran_number(0.0D0,1.0D0)
1074           xxh=betbol*(etot-etot0)
1075           if (xxh.lt.50.0D0) then
1076             xxh=dexp(-xxh)
1077             if (xxh.gt.xxr) accepted=.true.
1078           endif
1079        endif
1080        accepted=.true.
1081 !       print *,etot0,etot,accepted
1082        if (accepted) then
1083           etot0=etot
1084           call rmsd(rms)
1085           write(iout,*) 'etot=',i,etot,rms
1086           call write_pdb(i,'MC structure',etot)
1087 ! minimize
1088 !        call geom_to_var(nvar,var1)
1089         call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1090         call geom_to_var(nvar,var)
1091         call minimize(etot,var,iretcode,nfun)
1092         write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1093         call var_to_geom(nvar,var)
1094         call chainbuild
1095         call rmsd(rms)
1096         write(iout,*) 'etot mcm=',i,etot,rms
1097         call write_pdb(i+1,'MCM structure',etot)
1098         call var_to_geom(nvar,var1)
1099 ! --------
1100        else
1101           phi(jr)=phiold
1102        endif
1103       enddo
1104
1105 ! minimize
1106 !       call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1107 !       call geom_to_var(nvar,var)
1108 !
1109 !       call chainbuild
1110 !       call write_pdb(998 ,'sc min',etot)
1111 !
1112 !       call minimize(etot,var,iretcode,nfun)
1113 !       write(iout,*)'------------------------------------------------'
1114 !       write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
1115 !
1116 !       call var_to_geom(nvar,var)
1117 !       call chainbuild
1118 !       call write_pdb(999,'full min',etot)
1119
1120       return
1121       end subroutine test
1122 !-----------------------------------------------------------------------------
1123 !el#ifdef MPI
1124       subroutine test_n16
1125
1126 !el      use minim
1127 !      implicit real*8 (a-h,o-z)
1128 !      include 'DIMENSIONS'
1129       include 'mpif.h'
1130 !      include 'COMMON.GEO'
1131 !      include 'COMMON.VAR'
1132 !      include 'COMMON.INTERACT'
1133 !      include 'COMMON.IOUNITS'
1134 !      include 'COMMON.DISTFIT'
1135 !      include 'COMMON.SBRIDGE'
1136 !      include 'COMMON.CONTROL'
1137 !      include 'COMMON.FFIELD'
1138 !      include 'COMMON.MINIM'
1139 !      include 'COMMON.CHAIN'
1140       real(kind=8) :: time0,time1
1141       real(kind=8) :: energy(0:n_ene),ee
1142       real(kind=8),dimension(:),allocatable :: var,var1 !(maxvar) (maxvar=6*maxres)
1143       integer :: jdata(5)
1144       logical :: debug
1145 !el local variables
1146       integer :: i,ij,ieval,iretcode,nfun
1147       real(kind=8) :: etot
1148       debug=.true.
1149       allocate(var(6*nres),var1(6*nres))        !(maxvar) (maxvar=6*maxres)
1150 !
1151       call geom_to_var(nvar,var1)
1152       call chainbuild
1153       call etotal(energy)
1154       etot=energy(0)
1155       write(iout,*) nnt,nct,etot
1156       call write_pdb(1,'first structure',etot)
1157       call secondary2(.true.)
1158
1159       do i=1,4
1160         jdata(i)=bfrag(i,2)
1161       enddo
1162
1163       DO ij=1,4
1164        ieval=0
1165        jdata(5)=ij
1166        call var_to_geom(nvar,var1)
1167        write(iout,*) 'N16 test',(jdata(i),i=1,5)
1168        call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5), &
1169                       ieval,ij)
1170        call geom_to_var(nvar,var)
1171
1172       if (minim) then
1173 !el#ifdef MPI
1174        time0=MPI_WTIME()
1175 !el#endif
1176        call minimize(etot,var,iretcode,nfun)
1177        write(iout,*)'------------------------------------------------'
1178        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1179         '+ DIST eval',ieval
1180
1181 !el#ifdef MPI
1182        time1=MPI_WTIME()
1183 !el#endif
1184        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,&
1185                nfun/(time1-time0),' eval/s'
1186
1187        call var_to_geom(nvar,var)
1188        call chainbuild
1189        call write_pdb(ij*100+99,'full min',etot)
1190       endif
1191
1192
1193       ENDDO
1194
1195       return
1196       end subroutine test_n16
1197 !el#endif
1198 !-----------------------------------------------------------------------------
1199       subroutine test_local
1200
1201 !      implicit real*8 (a-h,o-z)
1202 !      include 'DIMENSIONS'
1203 !      include 'COMMON.GEO'
1204 !      include 'COMMON.VAR'
1205 !      include 'COMMON.INTERACT'
1206 !      include 'COMMON.IOUNITS'
1207       real(kind=8) :: time0,time1
1208       real(kind=8) :: energy(0:n_ene),ee
1209       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
1210       integer :: nft_sc
1211       real(kind=8) :: etot
1212 !
1213 !      allocate(varia(6*nres))  !(maxvar) (maxvar=6*maxres)
1214       call chainbuild
1215 !      call geom_to_var(nvar,varia)
1216       call write_pdb(1,'first structure',0d0)
1217
1218       call etotal(energy)
1219       etot=energy(0)
1220       write(iout,*) nnt,nct,etot
1221
1222       write(iout,*) 'calling sc_move'
1223       call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1224       write(iout,*) nft_sc,etot
1225       call write_pdb(2,'second structure',etot)
1226
1227       write(iout,*) 'calling local_move'
1228       call local_move_init(.false.)
1229       call local_move(24,29,20d0,50d0)
1230       call chainbuild
1231       call write_pdb(3,'third structure',etot)
1232
1233       write(iout,*) 'calling sc_move'
1234       call sc_move(24,29,5,10d0,nft_sc,etot)
1235       write(iout,*) nft_sc,etot
1236       call write_pdb(2,'last structure',etot)
1237
1238       return
1239       end subroutine test_local
1240 !-----------------------------------------------------------------------------
1241       subroutine test_sc
1242
1243 !      implicit real*8 (a-h,o-z)
1244 !      include 'DIMENSIONS'
1245 !      include 'COMMON.GEO'
1246 !      include 'COMMON.VAR'
1247 !      include 'COMMON.INTERACT'
1248 !      include 'COMMON.IOUNITS'
1249       real(kind=8) :: time0,time1,etot
1250       real(kind=8) :: energy(0:n_ene),ee
1251       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
1252       integer :: nft_sc
1253 !
1254       call chainbuild
1255 !      call geom_to_var(nvar,varia)
1256       call write_pdb(1,'first structure',0d0)
1257
1258       call etotal(energy)
1259       etot=energy(0)
1260       write(iout,*) nnt,nct,etot
1261
1262       write(iout,*) 'calling sc_move'
1263
1264       call sc_move(nnt,nct,5,10d0,nft_sc,etot)
1265       write(iout,*) nft_sc,etot
1266       call write_pdb(2,'second structure',etot)
1267
1268       write(iout,*) 'calling sc_move 2nd time'
1269
1270       call sc_move(nnt,nct,5,1d0,nft_sc,etot)
1271       write(iout,*) nft_sc,etot
1272       call write_pdb(3,'last structure',etot)
1273       return
1274       end subroutine test_sc
1275 !-----------------------------------------------------------------------------
1276       subroutine bgrow(bstrand,nbstrand,in,ind,new)
1277
1278 !      implicit real*8 (a-h,o-z)
1279 !      include 'DIMENSIONS'
1280 !      include 'COMMON.CHAIN'
1281       integer,dimension(nres/3,6) :: bstrand    !(maxres/3,6)
1282
1283 !el local variables
1284       integer :: nbstrand,in,ind,new,ishift,i
1285
1286       ishift=iabs(bstrand(in,ind+4)-new)
1287
1288       print *,'bgrow',bstrand(in,ind+4),new,ishift
1289
1290       bstrand(in,ind)=new
1291
1292       if(ind.eq.1)then
1293         bstrand(nbstrand,5)=bstrand(nbstrand,1)
1294         do i=1,nbstrand-1
1295           IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1296           if (bstrand(i,5).lt.bstrand(i,6)) then
1297             bstrand(i,5)=bstrand(i,5)-ishift
1298           else
1299             bstrand(i,5)=bstrand(i,5)+ishift
1300           endif
1301           ENDIF
1302         enddo
1303       else
1304         bstrand(nbstrand,6)=bstrand(nbstrand,2)
1305         do i=1,nbstrand-1
1306           IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
1307           if (bstrand(i,6).lt.bstrand(i,5)) then
1308             bstrand(i,6)=bstrand(i,6)-ishift
1309           else
1310             bstrand(i,6)=bstrand(i,6)+ishift
1311           endif
1312           ENDIF
1313         enddo
1314       endif
1315
1316       return
1317       end subroutine bgrow
1318 !-----------------------------------------------------------------------------
1319       subroutine test11
1320
1321       use geometry, only:dist
1322 !el      use minim
1323 !      implicit real*8 (a-h,o-z)
1324 !      include 'DIMENSIONS'
1325       include 'mpif.h'
1326 !      include 'COMMON.GEO'
1327 !      include 'COMMON.CHAIN'
1328 !      include 'COMMON.IOUNITS'
1329 !      include 'COMMON.VAR'
1330 !      include 'COMMON.CONTROL'
1331 !      include 'COMMON.SBRIDGE'
1332 !      include 'COMMON.FFIELD'
1333 !      include 'COMMON.MINIM'
1334 !
1335 !      include 'COMMON.DISTFIT'
1336       integer :: if(20,nres),nif,ifa(20)
1337       integer :: ibc(0:nres,0:nres),istrand(20)
1338       integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
1339       integer :: itmp(20,nres)
1340       real(kind=8) :: time0,time1
1341       real(kind=8) :: energy(0:n_ene),ee
1342       real(kind=8),dimension(6*nres) :: varia,vorg      !(maxvar) (maxvar=6*maxres)
1343 !
1344       logical :: debug,ltest,usedbfrag(nres/3)
1345       character(len=50) :: linia
1346 !
1347       integer :: betasheet(nres),ibetasheet(nres),nbetasheet
1348       integer :: bstrand(nres/3,6),nbstrand
1349       real(kind=8) :: etot
1350       integer :: i,j,jk,k,isa,m,l,ig,iconf,is,ii,iused_nbfrag,&
1351             in,ind,ifun,nfun,iretcode
1352 !------------------------
1353
1354       debug=.true.
1355 !------------------------
1356       nbstrand=0
1357       nbetasheet=0
1358       do i=1,nres
1359         betasheet(i)=0
1360         ibetasheet(i)=0
1361       enddo
1362       call geom_to_var(nvar,vorg)
1363       call secondary2(debug)
1364
1365       if (nbfrag.le.1) return
1366
1367       do i=1,nbfrag
1368          usedbfrag(i)=.false.
1369       enddo
1370
1371
1372       nbetasheet=nbetasheet+1
1373       nbstrand=2
1374       bstrand(1,1)=bfrag(1,1)
1375       bstrand(1,2)=bfrag(2,1)
1376       bstrand(1,3)=nbetasheet
1377       bstrand(1,4)=1
1378       bstrand(1,5)=bfrag(1,1)
1379       bstrand(1,6)=bfrag(2,1)
1380       do i=bfrag(1,1),bfrag(2,1)
1381         betasheet(i)=nbetasheet
1382         ibetasheet(i)=1
1383       enddo
1384 !
1385       bstrand(2,1)=bfrag(3,1)
1386       bstrand(2,2)=bfrag(4,1)
1387       bstrand(2,3)=nbetasheet
1388       bstrand(2,5)=bfrag(3,1)
1389       bstrand(2,6)=bfrag(4,1)
1390
1391       if (bfrag(3,1).le.bfrag(4,1)) then
1392         bstrand(2,4)=2
1393         do i=bfrag(3,1),bfrag(4,1)
1394           betasheet(i)=nbetasheet
1395           ibetasheet(i)=2
1396         enddo
1397       else
1398         bstrand(2,4)=-2
1399         do i=bfrag(4,1),bfrag(3,1)
1400           betasheet(i)=nbetasheet
1401           ibetasheet(i)=2
1402         enddo
1403       endif
1404
1405       iused_nbfrag=1
1406
1407       do while (iused_nbfrag.ne.nbfrag)
1408
1409       do j=2,nbfrag
1410
1411         IF (.not.usedbfrag(j)) THEN
1412
1413         write (*,*) j,(bfrag(i,j),i=1,4)
1414         do jk=6,1,-1
1415          write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand)
1416         enddo
1417         write (*,*) '------------------'
1418
1419
1420         if (bfrag(3,j).le.bfrag(4,j)) then
1421          do i=bfrag(3,j),bfrag(4,j)
1422           if(betasheet(i).eq.nbetasheet) then
1423             in=ibetasheet(i)
1424             do k=bfrag(3,j),bfrag(4,j)
1425               betasheet(k)=nbetasheet
1426               ibetasheet(k)=in
1427             enddo
1428             nbstrand=nbstrand+1
1429             usedbfrag(j)=.true.
1430             iused_nbfrag=iused_nbfrag+1
1431             do k=bfrag(1,j),bfrag(2,j)
1432               betasheet(k)=nbetasheet
1433               ibetasheet(k)=nbstrand
1434             enddo
1435             if (bstrand(in,4).lt.0) then
1436               bstrand(nbstrand,1)=bfrag(2,j)
1437               bstrand(nbstrand,2)=bfrag(1,j)
1438               bstrand(nbstrand,3)=nbetasheet
1439               bstrand(nbstrand,4)=-nbstrand
1440               bstrand(nbstrand,5)=bstrand(nbstrand,1)
1441               bstrand(nbstrand,6)=bstrand(nbstrand,2)
1442               if(bstrand(in,1).lt.bfrag(4,j)) then
1443                  call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1444               else
1445                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1446                      (bstrand(in,5)-bfrag(4,j))
1447               endif
1448               if(bstrand(in,2).gt.bfrag(3,j)) then
1449                  call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1450               else
1451                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1452                      (-bstrand(in,6)+bfrag(3,j))
1453               endif
1454             else
1455               bstrand(nbstrand,1)=bfrag(1,j)
1456               bstrand(nbstrand,2)=bfrag(2,j)
1457               bstrand(nbstrand,3)=nbetasheet
1458               bstrand(nbstrand,4)=nbstrand
1459               bstrand(nbstrand,5)=bstrand(nbstrand,1)
1460               bstrand(nbstrand,6)=bstrand(nbstrand,2)
1461               if(bstrand(in,1).gt.bfrag(3,j)) then
1462                  call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1463               else
1464                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1465                       (-bstrand(in,5)+bfrag(3,j))
1466               endif
1467               if(bstrand(in,2).lt.bfrag(4,j)) then
1468                  call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1469               else
1470                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1471                     (bstrand(in,6)-bfrag(4,j))
1472               endif
1473             endif
1474             goto 11
1475           endif
1476           if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then
1477             in=ibetasheet(bfrag(1,j)+i-bfrag(3,j))
1478             do k=bfrag(1,j),bfrag(2,j)
1479               betasheet(k)=nbetasheet
1480               ibetasheet(k)=in
1481             enddo
1482             nbstrand=nbstrand+1
1483             usedbfrag(j)=.true.
1484             iused_nbfrag=iused_nbfrag+1
1485             do k=bfrag(3,1),bfrag(4,1)
1486               betasheet(k)=nbetasheet
1487               ibetasheet(k)=nbstrand
1488             enddo
1489             if (bstrand(in,4).lt.0) then
1490               bstrand(nbstrand,1)=bfrag(4,j)
1491               bstrand(nbstrand,2)=bfrag(3,j)
1492               bstrand(nbstrand,3)=nbetasheet
1493               bstrand(nbstrand,4)=-nbstrand
1494               bstrand(nbstrand,5)=bstrand(nbstrand,1)
1495               bstrand(nbstrand,6)=bstrand(nbstrand,2)
1496               if(bstrand(in,1).lt.bfrag(2,j)) then
1497                  call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1498               else
1499                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1500                     (bstrand(in,5)-bfrag(2,j))
1501               endif
1502               if(bstrand(in,2).gt.bfrag(1,j)) then
1503                  call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1504               else
1505                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1506                     (-bstrand(in,6)+bfrag(1,j))
1507               endif
1508             else
1509               bstrand(nbstrand,1)=bfrag(3,j)
1510               bstrand(nbstrand,2)=bfrag(4,j)
1511               bstrand(nbstrand,3)=nbetasheet
1512               bstrand(nbstrand,4)=nbstrand
1513               bstrand(nbstrand,5)=bstrand(nbstrand,1)
1514               bstrand(nbstrand,6)=bstrand(nbstrand,2)
1515               if(bstrand(in,1).gt.bfrag(1,j)) then
1516                  call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1517               else
1518                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1519                     (-bstrand(in,5)+bfrag(1,j))
1520               endif
1521               if(bstrand(in,2).lt.bfrag(2,j)) then
1522                  call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1523               else
1524                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1525                      (bstrand(in,6)-bfrag(2,j))
1526               endif
1527             endif
1528             goto 11
1529           endif
1530          enddo
1531         else
1532          do i=bfrag(4,j),bfrag(3,j)
1533           if(betasheet(i).eq.nbetasheet) then
1534             in=ibetasheet(i)
1535             do k=bfrag(4,j),bfrag(3,j)
1536               betasheet(k)=nbetasheet
1537               ibetasheet(k)=in
1538             enddo
1539             nbstrand=nbstrand+1
1540             usedbfrag(j)=.true.
1541             iused_nbfrag=iused_nbfrag+1
1542             do k=bfrag(1,j),bfrag(2,j)
1543               betasheet(k)=nbetasheet
1544               ibetasheet(k)=nbstrand
1545             enddo
1546             if (bstrand(in,4).lt.0) then
1547               bstrand(nbstrand,1)=bfrag(1,j)
1548               bstrand(nbstrand,2)=bfrag(2,j)
1549               bstrand(nbstrand,3)=nbetasheet
1550               bstrand(nbstrand,4)=nbstrand
1551               bstrand(nbstrand,5)=bstrand(nbstrand,1)
1552               bstrand(nbstrand,6)=bstrand(nbstrand,2)
1553               if(bstrand(in,1).lt.bfrag(3,j)) then
1554                  call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
1555               else
1556                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1557                     (bstrand(in,5)-bfrag(3,j))
1558               endif
1559               if(bstrand(in,2).gt.bfrag(4,j)) then
1560                  call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
1561               else
1562                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1563                     (-bstrand(in,6)+bfrag(4,j))
1564               endif
1565             else
1566               bstrand(nbstrand,1)=bfrag(2,j)
1567               bstrand(nbstrand,2)=bfrag(1,j)
1568               bstrand(nbstrand,3)=nbetasheet
1569               bstrand(nbstrand,4)=-nbstrand
1570               bstrand(nbstrand,5)=bstrand(nbstrand,1)
1571               bstrand(nbstrand,6)=bstrand(nbstrand,2)
1572               if(bstrand(in,1).gt.bfrag(4,j)) then
1573                  call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
1574               else
1575                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1576                     (-bstrand(in,5)+bfrag(4,j))
1577               endif
1578               if(bstrand(in,2).lt.bfrag(3,j)) then
1579                  call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
1580               else
1581                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1582                    (bstrand(in,6)-bfrag(3,j))
1583               endif
1584             endif
1585             goto 11
1586           endif
1587           if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then
1588             in=ibetasheet(bfrag(2,j)-i+bfrag(4,j))
1589             do k=bfrag(1,j),bfrag(2,j)
1590               betasheet(k)=nbetasheet
1591               ibetasheet(k)=in
1592             enddo
1593             nbstrand=nbstrand+1
1594             usedbfrag(j)=.true.
1595             iused_nbfrag=iused_nbfrag+1
1596             do k=bfrag(4,j),bfrag(3,j)
1597               betasheet(k)=nbetasheet
1598               ibetasheet(k)=nbstrand
1599             enddo
1600             if (bstrand(in,4).lt.0) then
1601               bstrand(nbstrand,1)=bfrag(4,j)
1602               bstrand(nbstrand,2)=bfrag(3,j)
1603               bstrand(nbstrand,3)=nbetasheet
1604               bstrand(nbstrand,4)=nbstrand
1605               bstrand(nbstrand,5)=bstrand(nbstrand,1)
1606               bstrand(nbstrand,6)=bstrand(nbstrand,2)
1607               if(bstrand(in,1).lt.bfrag(2,j)) then
1608                  call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
1609               else
1610                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-&
1611                    (bstrand(in,5)-bfrag(2,j))
1612               endif
1613               if(bstrand(in,2).gt.bfrag(1,j)) then
1614                  call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
1615               else
1616                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+&
1617                    (-bstrand(in,6)+bfrag(1,j))
1618               endif
1619             else
1620               bstrand(nbstrand,1)=bfrag(3,j)
1621               bstrand(nbstrand,2)=bfrag(4,j)
1622               bstrand(nbstrand,3)=nbetasheet
1623               bstrand(nbstrand,4)=-nbstrand
1624               bstrand(nbstrand,5)=bstrand(nbstrand,1)
1625               bstrand(nbstrand,6)=bstrand(nbstrand,2)
1626               if(bstrand(in,1).gt.bfrag(1,j)) then
1627                  call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
1628               else
1629                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+&
1630                      (-bstrand(in,5)+bfrag(1,j))
1631               endif
1632               if(bstrand(in,2).lt.bfrag(2,j)) then
1633                  call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
1634               else
1635                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-&
1636                       (bstrand(in,6)-bfrag(2,j))
1637               endif
1638             endif
1639             goto 11
1640           endif
1641          enddo
1642         endif
1643
1644
1645
1646         ENDIF
1647       enddo
1648
1649       j=2
1650       do while (usedbfrag(j))
1651         j=j+1
1652       enddo
1653
1654       nbstrand=nbstrand+1
1655       nbetasheet=nbetasheet+1
1656       bstrand(nbstrand,1)=bfrag(1,j)
1657       bstrand(nbstrand,2)=bfrag(2,j)
1658       bstrand(nbstrand,3)=nbetasheet
1659       bstrand(nbstrand,5)=bfrag(1,j)
1660       bstrand(nbstrand,6)=bfrag(2,j)
1661
1662       bstrand(nbstrand,4)=nbstrand
1663       do i=bfrag(1,j),bfrag(2,j)
1664           betasheet(i)=nbetasheet
1665           ibetasheet(i)=nbstrand
1666       enddo
1667 !
1668       nbstrand=nbstrand+1
1669       bstrand(nbstrand,1)=bfrag(3,j)
1670       bstrand(nbstrand,2)=bfrag(4,j)
1671       bstrand(nbstrand,3)=nbetasheet
1672       bstrand(nbstrand,5)=bfrag(3,j)
1673       bstrand(nbstrand,6)=bfrag(4,j)
1674
1675       if (bfrag(3,j).le.bfrag(4,j)) then
1676         bstrand(nbstrand,4)=nbstrand
1677         do i=bfrag(3,j),bfrag(4,j)
1678           betasheet(i)=nbetasheet
1679           ibetasheet(i)=nbstrand
1680         enddo
1681       else
1682         bstrand(nbstrand,4)=-nbstrand
1683         do i=bfrag(4,j),bfrag(3,j)
1684           betasheet(i)=nbetasheet
1685           ibetasheet(i)=nbstrand
1686         enddo
1687       endif
1688
1689       iused_nbfrag=iused_nbfrag+1
1690       usedbfrag(j)=.true.
1691
1692
1693   11  continue
1694         do jk=6,1,-1
1695          write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand)
1696         enddo
1697
1698
1699       enddo
1700
1701       do i=1,nres
1702        if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i)
1703       enddo
1704       write(*,*)
1705       do j=6,1,-1
1706         write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand)
1707       enddo
1708
1709 !------------------------
1710       nifb=0
1711       do i=1,nbstrand
1712         do j=i+1,nbstrand
1713            if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or. &
1714                 iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then
1715              nifb=nifb+1
1716              ifb(nifb,1)=bstrand(i,4)
1717              ifb(nifb,2)=bstrand(j,4)
1718            endif
1719         enddo
1720       enddo
1721
1722       write(*,*)
1723       do i=1,nifb
1724          write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2)
1725       enddo
1726
1727       do i=1,nbstrand
1728            ifa(i)=bstrand(i,4)
1729       enddo
1730       write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand)
1731
1732       nif=iabs(bstrand(1,6)-bstrand(1,5))+1
1733       do j=2,nbstrand
1734        if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif) &
1735           nif=iabs(bstrand(j,6)-bstrand(j,5))+1
1736       enddo
1737
1738       write(*,*) nif
1739       do i=1,nif
1740         do j=1,nbstrand
1741           if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6))
1742           if (if(j,i).gt.0) then
1743             if(betasheet(if(j,i)).eq.0 .or. &
1744                 ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0
1745           else
1746             if(j,i)=0
1747           endif
1748         enddo
1749         write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand)
1750       enddo
1751
1752 !      read (inp,*) (ifa(i),i=1,4)
1753 !      do i=1,nres
1754 !       read (inp,*,err=20,end=20) (if(j,i),j=1,4)
1755 !      enddo
1756 ! 20   nif=i-1
1757        stop
1758 !------------------------
1759
1760       isa=4
1761       is=2*isa-1
1762       iconf=0
1763 !ccccccccccccccccccccccccccccccccc
1764       DO ig=1,is**isa-1
1765 !ccccccccccccccccccccccccccccccccc
1766
1767          ii=ig
1768          do j=1,is
1769            istrand(is-j+1)=int(ii/is**(is-j))
1770            ii=ii-istrand(is-j+1)*is**(is-j)
1771          enddo
1772          ltest=.true.
1773          do k=1,isa
1774            istrand(k)=istrand(k)+1
1775            if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1
1776          enddo
1777          do k=1,isa
1778            do l=1,isa
1779             if(istrand(k).eq.istrand(l).and.k.ne.l.or. &
1780                 istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false.
1781            enddo
1782          enddo
1783
1784          lifb0=1
1785          do m=1,nifb
1786            lifb(m)=0
1787            do k=1,isa-1
1788             if( &
1789                ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1790                ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1791              -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1792              -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1793                lifb(m)=1
1794            enddo
1795          lifb0=lifb0*lifb(m)
1796          enddo
1797
1798          if (mod(isa,2).eq.0) then
1799           do k=isa/2+1,isa
1800            if (istrand(k).eq.1) ltest=.false.
1801           enddo
1802          else
1803           do k=(isa+1)/2+1,isa
1804            if (istrand(k).eq.1) ltest=.false.
1805           enddo
1806          endif
1807
1808          IF (ltest.and.lifb0.eq.1) THEN
1809              iconf=iconf+1
1810
1811              call var_to_geom(nvar,vorg)
1812
1813              write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1814              write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
1815              write (linia,'(10i3)') (istrand(k),k=1,isa)
1816
1817       do i=1,nres
1818         do j=1,nres
1819          ibc(i,j)=0
1820         enddo
1821       enddo
1822
1823
1824       do i=1,4
1825        if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then
1826         do j=1,nif
1827          itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j)
1828         enddo
1829        else
1830         do j=1,nif
1831         itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1)
1832         enddo
1833        endif
1834       enddo
1835
1836       do i=1,nif
1837         write(*,*) (itmp(j,i),j=1,4)
1838       enddo
1839
1840       do i=1,nif
1841 !       ifa(1),ifa(2),ifa(3),ifa(4)
1842 !       if(1,i),if(2,i),if(3,i),if(4,i)
1843         do k=1,isa-1
1844          ltest=.false.
1845          do m=1,nifb
1846            if( &
1847                ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. &
1848                ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. &
1849              -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. &
1850              -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) &
1851          then
1852            ltest=.true.
1853            goto 110
1854          endif
1855          enddo
1856   110     continue
1857          if (ltest) then
1858           ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1
1859          else
1860           ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2
1861          endif
1862 !
1863         if (k.lt.3) &
1864          ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3
1865         if (k.lt.2) &
1866          ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4
1867         enddo
1868       enddo
1869 !------------------------
1870
1871 !
1872 !  freeze sec.elements
1873 !
1874        do i=1,nres
1875          mask(i)=1
1876          mask_phi(i)=1
1877          mask_theta(i)=1
1878          mask_side(i)=1
1879        enddo
1880
1881        do j=1,nbfrag
1882         do i=bfrag(1,j),bfrag(2,j)
1883          mask(i)=0
1884          mask_phi(i)=0
1885          mask_theta(i)=0
1886         enddo
1887         if (bfrag(3,j).le.bfrag(4,j)) then
1888          do i=bfrag(3,j),bfrag(4,j)
1889           mask(i)=0
1890           mask_phi(i)=0
1891           mask_theta(i)=0
1892          enddo
1893         else
1894          do i=bfrag(4,j),bfrag(3,j)
1895           mask(i)=0
1896           mask_phi(i)=0
1897           mask_theta(i)=0
1898          enddo
1899         endif
1900        enddo
1901        do j=1,nhfrag
1902         do i=hfrag(1,j),hfrag(2,j)
1903          mask(i)=0
1904          mask_phi(i)=0
1905          mask_theta(i)=0
1906         enddo
1907        enddo
1908        mask_r=.true.
1909
1910 !------------------------
1911 !      generate constrains
1912 !
1913        nhpb0=nhpb
1914        call chainbuild
1915        ind=0
1916        do i=1,nres-3
1917          do j=i+3,nres
1918           ind=ind+1
1919           if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then
1920             d0(ind)=DIST(i,j)
1921             w(ind)=10.0
1922             nhpb=nhpb+1
1923             ihpb(nhpb)=i
1924             jhpb(nhpb)=j
1925             forcon(nhpb)=10.0
1926             dhpb(nhpb)=d0(ind)
1927           else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then
1928             d0(ind)=5.0
1929             w(ind)=10.0
1930             nhpb=nhpb+1
1931             ihpb(nhpb)=i
1932             jhpb(nhpb)=j
1933             forcon(nhpb)=10.0
1934             dhpb(nhpb)=d0(ind)
1935           else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then
1936             d0(ind)=11.0
1937             w(ind)=10.0
1938             nhpb=nhpb+1
1939             ihpb(nhpb)=i
1940             jhpb(nhpb)=j
1941             forcon(nhpb)=10.0
1942             dhpb(nhpb)=d0(ind)
1943           else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then
1944             d0(ind)=16.0
1945             w(ind)=10.0
1946             nhpb=nhpb+1
1947             ihpb(nhpb)=i
1948             jhpb(nhpb)=j
1949             forcon(nhpb)=10.0
1950             dhpb(nhpb)=d0(ind)
1951           else if ( ibc(i,j).gt.0 ) then
1952             d0(ind)=DIST(i,ibc(i,j))
1953             w(ind)=10.0
1954             nhpb=nhpb+1
1955             ihpb(nhpb)=i
1956             jhpb(nhpb)=j
1957             forcon(nhpb)=10.0
1958             dhpb(nhpb)=d0(ind)
1959           else if ( ibc(j,i).gt.0 ) then
1960             d0(ind)=DIST(ibc(j,i),j)
1961             w(ind)=10.0
1962             nhpb=nhpb+1
1963             ihpb(nhpb)=i
1964             jhpb(nhpb)=j
1965             forcon(nhpb)=10.0
1966             dhpb(nhpb)=d0(ind)
1967           else
1968             w(ind)=0.0
1969           endif
1970           ddd(ind)=d0(ind)
1971          enddo
1972        enddo
1973        call hpb_partition
1974 !d--------------------------
1975
1976       write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
1977               ibc(jhpb(i),ihpb(i)),' --',&
1978               ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
1979
1980 !d      nhpb=0
1981 !d      goto 901
1982 !
1983 !
1984 !el#ifdef MPI
1985       call contact_cp_min(varia,ifun,iconf,linia,debug)
1986       if (minim) then
1987 !el#ifdef MPI
1988        time0=MPI_WTIME()
1989 !el#endif
1990        call minimize(etot,varia,iretcode,nfun)
1991        write(iout,*)'------------------------------------------------'
1992        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
1993         '+ DIST eval',ifun
1994
1995 !el#ifdef MPI
1996        time1=MPI_WTIME()
1997 !el#endif
1998        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,&
1999                nfun/(time1-time0),' eval/s'
2000
2001        write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa)
2002        call var_to_geom(nvar,varia)
2003        call chainbuild
2004        call write_pdb(900+iconf,linia,etot)
2005       endif
2006 !el#endif
2007       call etotal(energy)
2008       etot=energy(0)
2009       call enerprint(energy)
2010 !d      call intout
2011 !d      call briefout(0,etot)
2012 !d      call secondary2(.true.)
2013
2014  901  CONTINUE
2015 !test      return
2016 !ccccccccccccccccccccccccccccccccccc
2017       ENDIF
2018       ENDDO
2019 !ccccccccccccccccccccccccccccccccccc
2020
2021       return
2022   10  write (iout,'(a)') 'Error reading test structure.'
2023       return
2024       end subroutine test11
2025 !-----------------------------------------------------------------------------
2026       subroutine test3
2027
2028       use geometry, only:dist
2029 !el      use minim
2030 !      implicit real*8 (a-h,o-z)
2031 !      include 'DIMENSIONS'
2032       include 'mpif.h'
2033 !      include 'COMMON.GEO'
2034 !      include 'COMMON.CHAIN'
2035 !      include 'COMMON.IOUNITS'
2036 !      include 'COMMON.VAR'
2037 !      include 'COMMON.CONTROL'
2038 !      include 'COMMON.SBRIDGE'
2039 !      include 'COMMON.FFIELD'
2040 !      include 'COMMON.MINIM'
2041 !
2042 !      include 'COMMON.DISTFIT'
2043       integer :: if(3,nres),nif
2044       integer :: ibc(nres,nres),istrand(20)
2045       integer :: ibd(nres),ifb(10,2),nifb,lifb(10),lifb0
2046       real(kind=8) :: time0,time1
2047       real(kind=8) :: energy(0:n_ene),ee
2048       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
2049 !
2050       logical :: debug,ltest
2051       character(len=50) :: linia
2052       integer :: ieval,i,j,ind,in_pdb,nfun,iretcode
2053       real(kind=8) :: etot
2054 !
2055       do i=1,nres
2056        read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i)
2057       enddo
2058  20   nif=i-1
2059       write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),&
2060                                                i=1,nif)
2061
2062
2063 !------------------------
2064       call secondary2(debug)
2065 !------------------------
2066       do i=1,nres
2067         do j=1,nres
2068          ibc(i,j)=0
2069         enddo
2070       enddo
2071
2072 !
2073 !  freeze sec.elements and store indexes for beta constrains
2074 !
2075        do i=1,nres
2076          mask(i)=1
2077          mask_phi(i)=1
2078          mask_theta(i)=1
2079          mask_side(i)=1
2080        enddo
2081
2082        do j=1,nbfrag
2083         do i=bfrag(1,j),bfrag(2,j)
2084          mask(i)=0
2085          mask_phi(i)=0
2086          mask_theta(i)=0
2087         enddo
2088         if (bfrag(3,j).le.bfrag(4,j)) then
2089          do i=bfrag(3,j),bfrag(4,j)
2090           mask(i)=0
2091           mask_phi(i)=0
2092           mask_theta(i)=0
2093           ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1
2094          enddo
2095         else
2096          do i=bfrag(4,j),bfrag(3,j)
2097           mask(i)=0
2098           mask_phi(i)=0
2099           mask_theta(i)=0
2100           ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1
2101          enddo
2102         endif
2103        enddo
2104        do j=1,nhfrag
2105         do i=hfrag(1,j),hfrag(2,j)
2106          mask(i)=0
2107          mask_phi(i)=0
2108          mask_theta(i)=0
2109         enddo
2110        enddo
2111        mask_r=.true.
2112
2113
2114 ! ---------------- test --------------
2115        do i=1,nif
2116          if (ibc(if(1,i),if(2,i)).eq.-1) then
2117              ibc(if(1,i),if(2,i))=if(3,i)
2118              ibc(if(1,i),if(3,i))=if(2,i)
2119          else if (ibc(if(2,i),if(1,i)).eq.-1) then
2120              ibc(if(2,i),if(1,i))=0
2121              ibc(if(1,i),if(2,i))=if(3,i)
2122              ibc(if(1,i),if(3,i))=if(2,i)
2123          else
2124              ibc(if(1,i),if(2,i))=if(3,i)
2125              ibc(if(1,i),if(3,i))=if(2,i)
2126          endif
2127        enddo
2128
2129        do i=1,nres
2130         do j=1,nres
2131          if (ibc(i,j).ne.0)  write(*,'(3i5)') i,j,ibc(i,j)
2132         enddo
2133        enddo
2134 !------------------------
2135        call chainbuild
2136        ind=0
2137        do i=1,nres-3
2138          do j=i+3,nres
2139           ind=ind+1
2140           if ( ibc(i,j).eq.-1 ) then
2141             d0(ind)=DIST(i,j)
2142             w(ind)=10.0
2143             nhpb=nhpb+1
2144             ihpb(nhpb)=i
2145             jhpb(nhpb)=j
2146             forcon(nhpb)=10.0
2147             dhpb(nhpb)=d0(ind)
2148           else if ( ibc(i,j).gt.0 ) then
2149             d0(ind)=DIST(i,ibc(i,j))
2150             w(ind)=10.0
2151             nhpb=nhpb+1
2152             ihpb(nhpb)=i
2153             jhpb(nhpb)=j
2154             forcon(nhpb)=10.0
2155             dhpb(nhpb)=d0(ind)
2156           else if ( ibc(j,i).gt.0 ) then
2157             d0(ind)=DIST(ibc(j,i),j)
2158             w(ind)=10.0
2159             nhpb=nhpb+1
2160             ihpb(nhpb)=i
2161             jhpb(nhpb)=j
2162             forcon(nhpb)=10.0
2163             dhpb(nhpb)=d0(ind)
2164           else
2165             w(ind)=0.0
2166           endif
2167          enddo
2168        enddo
2169        call hpb_partition
2170
2171 !d--------------------------
2172       write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),&
2173               ibc(jhpb(i),ihpb(i)),' --',&
2174               ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
2175
2176
2177       linia='dist'
2178       debug=.true.
2179       in_pdb=7
2180 !
2181 !el#ifdef MPI
2182       call contact_cp_min(varia,ieval,in_pdb,linia,debug)
2183       if (minim) then
2184 !el#ifdef MPI
2185        time0=MPI_WTIME()
2186 !el#endif
2187        call minimize(etot,varia,iretcode,nfun)
2188        write(iout,*)'------------------------------------------------'
2189        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2190         '+ DIST eval',ieval
2191
2192 !el#ifdef MPI
2193        time1=MPI_WTIME()
2194 !el#endif
2195        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,&
2196                nfun/(time1-time0),' eval/s'
2197
2198
2199        call var_to_geom(nvar,varia)
2200        call chainbuild
2201        call write_pdb(999,'full min',etot)
2202       endif
2203 !el#endif
2204       call etotal(energy)
2205       etot=energy(0)
2206       call enerprint(energy)
2207       call intout
2208       call briefout(0,etot)
2209       call secondary2(.true.)
2210
2211       return
2212   10  write (iout,'(a)') 'Error reading test structure.'
2213       return
2214       end subroutine test3
2215 !-----------------------------------------------------------------------------
2216       subroutine test__
2217
2218 !el      use minim
2219 !      implicit real*8 (a-h,o-z)
2220 !      include 'DIMENSIONS'
2221       include 'mpif.h'
2222 !      include 'COMMON.GEO'
2223 !      include 'COMMON.CHAIN'
2224 !      include 'COMMON.IOUNITS'
2225 !      include 'COMMON.VAR'
2226 !      include 'COMMON.CONTROL'
2227 !      include 'COMMON.SBRIDGE'
2228 !      include 'COMMON.FFIELD'
2229 !      include 'COMMON.MINIM'
2230 !
2231 !      include 'COMMON.DISTFIT'
2232       integer :: if(2,2),ind
2233       integer :: iff(nres)
2234       real(kind=8) :: time0,time1
2235       real(kind=8) :: energy(0:n_ene),ee
2236       real(kind=8),dimension(nres) :: theta2,phi2,alph2,omeg2,&
2237                                       theta1,phi1,alph1,omeg1   !(maxres)
2238       real(kind=8),dimension(6*nres) :: varia,varia2    !(maxvar) (maxvar=6*maxres)
2239 !
2240       integer :: i,j,nn,ifun,iretcode,nfun
2241       real(kind=8) :: etot
2242       nn=0
2243
2244       read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2)
2245       write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2)
2246       read (inp,*,err=10,end=10) (theta2(i),i=3,nres)
2247       read (inp,*,err=10,end=10) (phi2(i),i=4,nres)
2248       read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1)
2249       read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1)
2250       do i=1,nres
2251         theta2(i)=deg2rad*theta2(i)
2252         phi2(i)=deg2rad*phi2(i)
2253         alph2(i)=deg2rad*alph2(i)
2254         omeg2(i)=deg2rad*omeg2(i)
2255       enddo
2256       do i=1,nres
2257         theta1(i)=theta(i)
2258         phi1(i)=phi(i)
2259         alph1(i)=alph(i)
2260         omeg1(i)=omeg(i)
2261       enddo
2262
2263       do i=1,nres
2264        mask(i)=1
2265       enddo
2266
2267
2268 !------------------------
2269       do i=1,nres
2270          iff(i)=0
2271       enddo
2272       do j=1,2
2273         do i=if(j,1),if(j,2)
2274           iff(i)=1
2275         enddo
2276       enddo
2277
2278       call chainbuild
2279       call geom_to_var(nvar,varia)
2280       call write_pdb(1,'first structure',0d0)
2281
2282         call secondary(.true.)
2283
2284         call secondary2(.true.)
2285
2286         do j=1,nbfrag
2287            if ( (bfrag(3,j).lt.bfrag(4,j) .or. &
2288                 bfrag(4,j)-bfrag(2,j).gt.4) .and. &
2289                 bfrag(2,j)-bfrag(1,j).gt.3 ) then
2290              nn=nn+1
2291
2292              if (bfrag(3,j).lt.bfrag(4,j)) then
2293                write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2294                            "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2295                                 ",",bfrag(3,j)-1,"-",bfrag(4,j)-1
2296              else
2297                write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2298                            "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
2299                                 ",",bfrag(4,j)-1,"-",bfrag(3,j)-1
2300              endif
2301            endif
2302         enddo
2303
2304       do i=1,nres
2305         theta(i)=theta2(i)
2306         phi(i)=phi2(i)
2307         alph(i)=alph2(i)
2308         omeg(i)=omeg2(i)
2309       enddo
2310
2311       call chainbuild
2312       call geom_to_var(nvar,varia2)
2313       call write_pdb(2,'second structure',0d0)
2314
2315
2316
2317 !-------------------------------------------------------
2318 !el#ifdef MPI
2319       ifun=-1
2320       call contact_cp(varia,varia2,iff,ifun,7)
2321       if (minim) then
2322 !el#ifdef MPI
2323        time0=MPI_WTIME()
2324 !el#endif
2325        call minimize(etot,varia,iretcode,nfun)
2326        write(iout,*)'------------------------------------------------'
2327        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
2328         '+ DIST eval',ifun
2329
2330 !el#ifdef MPI
2331        time1=MPI_WTIME()
2332 !el#endif
2333        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,&
2334                nfun/(time1-time0),' eval/s'
2335
2336
2337        call var_to_geom(nvar,varia)
2338        call chainbuild
2339        call write_pdb(999,'full min',etot)
2340       endif
2341 !el#endif
2342       call etotal(energy)
2343       etot=energy(0)
2344       call enerprint(energy)
2345       call intout
2346       call briefout(0,etot)
2347
2348       return
2349   10  write (iout,'(a)') 'Error reading test structure.'
2350       return
2351       end subroutine test__
2352 !-----------------------------------------------------------------------------
2353       subroutine secondary(lprint)
2354
2355 !      implicit real*8 (a-h,o-z)
2356 !      include 'DIMENSIONS'
2357 !      include 'COMMON.CHAIN'
2358 !      include 'COMMON.IOUNITS'
2359 !      include 'COMMON.DISTFIT'
2360
2361       integer :: ncont,icont(2,nres*nres/2),isec(nres,3)
2362       logical :: lprint,not_done
2363       real(kind=4) :: dcont(nres*nres/2),d
2364       real(kind=4) :: rcomp = 7.0
2365       real(kind=4) :: rbeta = 5.2
2366       real(kind=4) :: ralfa = 5.2
2367       real(kind=4) :: r310 = 6.6
2368       real(kind=8),dimension(3) :: xpi,xpj
2369       integer :: i,k,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,iii1,jjj1,&
2370             nhelix
2371       call chainbuild_cart
2372 !d      call write_pdb(99,'sec structure',0d0)
2373       ncont=0
2374       nbfrag=0
2375       nhfrag=0
2376       do i=1,nres
2377         isec(i,1)=0
2378         isec(i,2)=0
2379         isec(i,3)=0
2380       enddo
2381
2382       do i=2,nres-3
2383         do k=1,3
2384           xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
2385         enddo
2386         do j=i+2,nres
2387           do k=1,3
2388              xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
2389           enddo
2390 !d        d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
2391 !d     &         (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
2392 !d     &         (c(3,i)-c(3,j))*(c(3,i)-c(3,j))
2393 !d          print *,'CA',i,j,d
2394           d =  (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) + &
2395                (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) + &
2396                (xpi(3)-xpj(3))*(xpi(3)-xpj(3))
2397           if ( d.lt.rcomp*rcomp) then
2398             ncont=ncont+1
2399             icont(1,ncont)=i
2400             icont(2,ncont)=j
2401             dcont(ncont)=sqrt(d)
2402           endif
2403         enddo
2404       enddo
2405       if (lprint) then
2406         write (iout,*)
2407         write (iout,'(a)') '#PP contact map distances:'
2408         do i=1,ncont
2409           write (iout,'(3i4,f10.5)') &
2410            i,icont(1,i),icont(2,i),dcont(i)
2411         enddo
2412       endif
2413
2414 ! finding parallel beta
2415 !d      write (iout,*) '------- looking for parallel beta -----------'
2416       nbeta=0
2417       nstrand=0
2418       do i=1,ncont
2419         i1=icont(1,i)
2420         j1=icont(2,i)
2421         if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and. &
2422             isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2423           (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2424           (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2425           (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2426           (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2427            ) then
2428           ii1=i1
2429           jj1=j1
2430 !d         write (iout,*) i1,j1,dcont(i)
2431           not_done=.true.
2432           do while (not_done)
2433            i1=i1+1
2434            j1=j1+1
2435             do j=1,ncont
2436               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) &
2437                     .and. dcont(j).le.rbeta .and. &
2438             isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2439           (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2440           (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2441           (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2442           (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2443                                   ) goto 5
2444             enddo
2445             not_done=.false.
2446   5         continue
2447 !d            write (iout,*) i1,j1,dcont(j),not_done
2448           enddo
2449           j1=j1-1
2450           i1=i1-1
2451           if (i1-ii1.gt.1) then
2452             ii1=max0(ii1-1,1)
2453             jj1=max0(jj1-1,1)
2454             nbeta=nbeta+1
2455             if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
2456
2457             nbfrag=nbfrag+1
2458             bfrag(1,nbfrag)=ii1
2459             bfrag(2,nbfrag)=i1
2460             bfrag(3,nbfrag)=jj1
2461             bfrag(4,nbfrag)=j1
2462
2463             do ij=ii1,i1
2464              isec(ij,1)=isec(ij,1)+1
2465              isec(ij,1+isec(ij,1))=nbeta
2466             enddo
2467             do ij=jj1,j1
2468              isec(ij,1)=isec(ij,1)+1
2469              isec(ij,1+isec(ij,1))=nbeta
2470             enddo
2471
2472            if(lprint) then
2473             nstrand=nstrand+1
2474             if (nbeta.le.9) then
2475               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2476                 "DefPropRes 'strand",nstrand,&
2477                 "' 'num = ",ii1-1,"..",i1-1,"'"
2478             else
2479               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2480                 "DefPropRes 'strand",nstrand,&
2481                 "' 'num = ",ii1-1,"..",i1-1,"'"
2482             endif
2483             nstrand=nstrand+1
2484             if (nbeta.le.9) then
2485               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2486                 "DefPropRes 'strand",nstrand,&
2487                 "' 'num = ",jj1-1,"..",j1-1,"'"
2488             else
2489               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2490                 "DefPropRes 'strand",nstrand,&
2491                 "' 'num = ",jj1-1,"..",j1-1,"'"
2492             endif
2493               write(12,'(a8,4i4)') &
2494                 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
2495            endif
2496           endif
2497         endif
2498       enddo
2499
2500 ! finding antiparallel beta
2501 !d      write (iout,*) '--------- looking for antiparallel beta ---------'
2502
2503       do i=1,ncont
2504         i1=icont(1,i)
2505         j1=icont(2,i)
2506         if (dcont(i).le.rbeta.and. &
2507             isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2508           (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2509           (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2510           (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2511           (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2512            ) then
2513           ii1=i1
2514           jj1=j1
2515 !d          write (iout,*) i1,j1,dcont(i)
2516
2517           not_done=.true.
2518           do while (not_done)
2519            i1=i1+1
2520            j1=j1-1
2521             do j=1,ncont
2522               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
2523             isec(i1,1).le.1.and.isec(j1,1).le.1.and. &
2524           (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. &
2525           (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. &
2526           (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. &
2527           (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0) &
2528                  .and. dcont(j).le.rbeta ) goto 6
2529             enddo
2530             not_done=.false.
2531   6         continue
2532 !d            write (iout,*) i1,j1,dcont(j),not_done
2533           enddo
2534           i1=i1-1
2535           j1=j1+1
2536           if (i1-ii1.gt.1) then
2537             if(lprint)write (iout,*)'antiparallel beta',&
2538                          nbeta,ii1-1,i1,jj1,j1-1
2539
2540             nbfrag=nbfrag+1
2541             bfrag(1,nbfrag)=max0(ii1-1,1)
2542             bfrag(2,nbfrag)=i1
2543             bfrag(3,nbfrag)=jj1
2544             bfrag(4,nbfrag)=max0(j1-1,1)
2545
2546             nbeta=nbeta+1
2547             iii1=max0(ii1-1,1)
2548             do ij=iii1,i1
2549              isec(ij,1)=isec(ij,1)+1
2550              isec(ij,1+isec(ij,1))=nbeta
2551             enddo
2552             jjj1=max0(j1-1,1)
2553             do ij=jjj1,jj1
2554              isec(ij,1)=isec(ij,1)+1
2555              isec(ij,1+isec(ij,1))=nbeta
2556             enddo
2557
2558
2559            if (lprint) then
2560             nstrand=nstrand+1
2561             if (nstrand.le.9) then
2562               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2563                 "DefPropRes 'strand",nstrand,&
2564                 "' 'num = ",ii1-2,"..",i1-1,"'"
2565             else
2566               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2567                 "DefPropRes 'strand",nstrand,&
2568                 "' 'num = ",ii1-2,"..",i1-1,"'"
2569             endif
2570             nstrand=nstrand+1
2571             if (nstrand.le.9) then
2572               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
2573                 "DefPropRes 'strand",nstrand,&
2574                 "' 'num = ",j1-2,"..",jj1-1,"'"
2575             else
2576               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
2577                 "DefPropRes 'strand",nstrand,&
2578                 "' 'num = ",j1-2,"..",jj1-1,"'"
2579             endif
2580               write(12,'(a8,4i4)') &
2581                 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
2582            endif
2583           endif
2584         endif
2585       enddo
2586
2587       if (nstrand.gt.0.and.lprint) then
2588         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
2589         do i=2,nstrand
2590          if (i.le.9) then
2591           write(12,'(a9,i1,$)') " | strand",i
2592          else
2593           write(12,'(a9,i2,$)') " | strand",i
2594          endif
2595         enddo
2596         write(12,'(a1)') "'"
2597       endif
2598
2599
2600 ! finding alpha or 310 helix
2601
2602       nhelix=0
2603       do i=1,ncont
2604         i1=icont(1,i)
2605         j1=icont(2,i)
2606         if (j1.eq.i1+3.and.dcont(i).le.r310 &
2607            .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
2608 !d          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
2609 !d          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
2610           ii1=i1
2611           jj1=j1
2612           if (isec(ii1,1).eq.0) then
2613             not_done=.true.
2614           else
2615             not_done=.false.
2616           endif
2617           do while (not_done)
2618             i1=i1+1
2619             j1=j1+1
2620             do j=1,ncont
2621               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
2622             enddo
2623             not_done=.false.
2624   10        continue
2625 !d            write (iout,*) i1,j1,not_done
2626           enddo
2627           j1=j1-1
2628           if (j1-ii1.gt.4) then
2629             nhelix=nhelix+1
2630 !d            write (iout,*)'helix',nhelix,ii1,j1
2631
2632             nhfrag=nhfrag+1
2633             hfrag(1,nhfrag)=ii1
2634             hfrag(2,nhfrag)=max0(j1-1,1)
2635
2636             do ij=ii1,j1
2637              isec(ij,1)=-1
2638             enddo
2639            if (lprint) then
2640             write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
2641             if (nhelix.le.9) then
2642               write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
2643                 "DefPropRes 'helix",nhelix,&
2644                 "' 'num = ",ii1-1,"..",j1-2,"'"
2645             else
2646               write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
2647                 "DefPropRes 'helix",nhelix,&
2648                 "' 'num = ",ii1-1,"..",j1-2,"'"
2649             endif
2650            endif
2651           endif
2652         endif
2653       enddo
2654
2655       if (nhelix.gt.0.and.lprint) then
2656         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
2657         do i=2,nhelix
2658          if (nhelix.le.9) then
2659           write(12,'(a8,i1,$)') " | helix",i
2660          else
2661           write(12,'(a8,i2,$)') " | helix",i
2662          endif
2663         enddo
2664         write(12,'(a1)') "'"
2665       endif
2666
2667       if (lprint) then
2668        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
2669        write(12,'(a20)') "XMacStand ribbon.mac"
2670       endif
2671
2672       return
2673       end subroutine secondary
2674 !-----------------------------------------------------------------------------
2675       subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
2676
2677       use geometry, only:dist
2678 !el      use minim
2679 !      implicit real*8 (a-h,o-z)
2680 !      include 'DIMENSIONS'
2681       include 'mpif.h'
2682 !      include 'COMMON.SBRIDGE'
2683 !      include 'COMMON.FFIELD'
2684 !      include 'COMMON.IOUNITS'
2685 !      include 'COMMON.DISTFIT'
2686 !      include 'COMMON.VAR'
2687 !      include 'COMMON.CHAIN'
2688 !      include 'COMMON.MINIM'
2689
2690       character(len=50) :: linia
2691       integer :: nf,ij(4)
2692       real(kind=8),dimension(6*nres) :: var,var2        !(maxvar) (maxvar=6*maxres)
2693       real(kind=8) :: time0,time1
2694       integer :: iff(nres),ieval
2695       real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1   !(maxres)
2696
2697 !el local variables
2698       integer :: in_pdb,i,j,ind,ipot0,maxmin0,maxfun0,nfun,iwsk,iretcode
2699       real(kind=8) :: wstrain0,etot
2700       integer :: maxres22
2701       maxres22=nres*(nres+1)/2
2702
2703       if(.not.allocated(DRDG)) allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
2704        call var_to_geom(nvar,var)
2705        call chainbuild
2706        nhpb0=nhpb
2707        ind=0
2708        do i=1,nres-3
2709          do j=i+3,nres
2710           ind=ind+1
2711           if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2712             d0(ind)=DIST(i,j)
2713             w(ind)=10.0
2714             nhpb=nhpb+1
2715             ihpb(nhpb)=i
2716             jhpb(nhpb)=j
2717             forcon(nhpb)=10.0
2718             dhpb(nhpb)=d0(ind)
2719           else
2720             w(ind)=0.0
2721           endif
2722          enddo
2723        enddo
2724        call hpb_partition
2725
2726        do i=1,nres
2727         theta1(i)=theta(i)
2728         phi1(i)=phi(i)
2729         alph1(i)=alph(i)
2730         omeg1(i)=omeg(i)
2731        enddo
2732
2733        call var_to_geom(nvar,var2)
2734
2735        do i=1,nres
2736           if ( iff(i).eq.1 ) then
2737               theta(i)=theta1(i)
2738               phi(i)=phi1(i)
2739               alph(i)=alph1(i)
2740               omeg(i)=omeg1(i)
2741           endif
2742        enddo
2743
2744        call chainbuild
2745 !d       call write_pdb(3,'combined structure',0d0)
2746 !d       time0=MPI_WTIME()
2747
2748        NX=NRES-3
2749        NY=((NRES-4)*(NRES-5))/2
2750        call distfit(.true.,200)
2751
2752 !d       time1=MPI_WTIME()
2753 !d       write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
2754
2755        ipot0=ipot
2756        maxmin0=maxmin
2757        maxfun0=maxfun
2758        wstrain0=wstrain
2759
2760        ipot=6
2761        maxmin=2000
2762        maxfun=5000
2763        call geom_to_var(nvar,var)
2764 !d       time0=MPI_WTIME()
2765        call minimize(etot,var,iretcode,nfun)
2766        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2767
2768 !d       time1=MPI_WTIME()
2769 !d       write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2770 !d     &         nfun/(time1-time0),' SOFT eval/s'
2771         call var_to_geom(nvar,var)
2772         call chainbuild
2773
2774
2775         iwsk=0
2776         nf=0
2777         if (iff(1).eq.1) then
2778           iwsk=1
2779           nf=nf+1
2780           ij(nf)=0
2781         endif
2782         do i=2,nres
2783            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2784              iwsk=1
2785              nf=nf+1
2786              ij(nf)=i
2787            endif
2788            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2789              iwsk=0
2790              nf=nf+1
2791              ij(nf)=i-1
2792            endif
2793         enddo
2794         if (iff(nres).eq.1) then
2795           nf=nf+1
2796           ij(nf)=nres
2797         endif
2798
2799
2800 !d        write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
2801 !d     &                     "select",ij(1),"-",ij(2),
2802 !d     &                         ",",ij(3),"-",ij(4)
2803 !d        call write_pdb(in_pdb,linia,etot)
2804
2805
2806        ipot=ipot0
2807        maxmin=maxmin0
2808        maxfun=maxfun0
2809 !d       time0=MPI_WTIME()
2810        call minimize(etot,var,iretcode,nfun)
2811 !d       write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2812        ieval=nfun
2813
2814 !d       time1=MPI_WTIME()
2815 !d       write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
2816 !d     &         nfun/(time1-time0),' eval/s'
2817 !d       call var_to_geom(nvar,var)
2818 !d       call chainbuild
2819 !d       call write_pdb(6,'dist structure',etot)
2820
2821
2822        nhpb= nhpb0
2823        link_start=1
2824        link_end=nhpb
2825        wstrain=wstrain0
2826
2827       return
2828       end subroutine contact_cp2
2829 !-----------------------------------------------------------------------------
2830       subroutine contact_cp(var,var2,iff,ieval,in_pdb)
2831
2832       use geometry, only:dist
2833 !el      use minim
2834 !      implicit real*8 (a-h,o-z)
2835 !      include 'DIMENSIONS'
2836 !      include 'COMMON.SBRIDGE'
2837 !      include 'COMMON.FFIELD'
2838 !      include 'COMMON.IOUNITS'
2839 !      include 'COMMON.DISTFIT'
2840 !      include 'COMMON.VAR'
2841 !      include 'COMMON.CHAIN'
2842 !      include 'COMMON.MINIM'
2843
2844       character(len=50) :: linia
2845       integer :: nf,ij(4)
2846       real(kind=8) :: energy(0:n_ene)
2847       real(kind=8),dimension(6*nres) :: var,var2        !(maxvar) (maxvar=6*maxres)
2848       real(kind=8) :: time0,time1
2849       integer :: iff(nres),ieval
2850       real(kind=8),dimension(nres) :: theta1,phi1,alph1,omeg1   !(maxres)
2851       logical :: debug
2852
2853 !el local variables
2854       integer :: in_pdb,i,j,ind,iwsk
2855
2856       debug=.false.
2857 !      debug=.true.
2858       if (ieval.eq.-1) debug=.true.
2859
2860
2861 !
2862 ! store selected dist. constrains from 1st structure
2863 !
2864 #ifdef OSF
2865 !     Intercept NaNs in the coordinates
2866 !      write(iout,*) (var(i),i=1,nvar)
2867       x_sum=0.D0
2868       do i=1,nvar
2869         x_sum=x_sum+var(i)
2870       enddo
2871       if (x_sum.ne.x_sum) then
2872         write(iout,*)" *** contact_cp : Found NaN in coordinates"
2873         call flush(iout)
2874         print *," *** contact_cp : Found NaN in coordinates"
2875         return
2876       endif
2877 #endif
2878
2879
2880        call var_to_geom(nvar,var)
2881        call chainbuild
2882        nhpb0=nhpb
2883        ind=0
2884        do i=1,nres-3
2885          do j=i+3,nres
2886           ind=ind+1
2887           if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2888             d0(ind)=DIST(i,j)
2889             w(ind)=10.0
2890             nhpb=nhpb+1
2891             ihpb(nhpb)=i
2892             jhpb(nhpb)=j
2893             forcon(nhpb)=10.0
2894             dhpb(nhpb)=d0(ind)
2895           else
2896             w(ind)=0.0
2897           endif
2898          enddo
2899        enddo
2900        call hpb_partition
2901
2902        do i=1,nres
2903         theta1(i)=theta(i)
2904         phi1(i)=phi(i)
2905         alph1(i)=alph(i)
2906         omeg1(i)=omeg(i)
2907        enddo
2908
2909 !
2910 !  freeze sec.elements from 2nd structure
2911 !
2912        do i=1,nres
2913          mask_phi(i)=1
2914          mask_theta(i)=1
2915          mask_side(i)=1
2916        enddo
2917
2918        call var_to_geom(nvar,var2)
2919        call secondary2(debug)
2920        do j=1,nbfrag
2921         do i=bfrag(1,j),bfrag(2,j)
2922          mask(i)=0
2923          mask_phi(i)=0
2924          mask_theta(i)=0
2925         enddo
2926         if (bfrag(3,j).le.bfrag(4,j)) then
2927          do i=bfrag(3,j),bfrag(4,j)
2928           mask(i)=0
2929           mask_phi(i)=0
2930           mask_theta(i)=0
2931          enddo
2932         else
2933          do i=bfrag(4,j),bfrag(3,j)
2934           mask(i)=0
2935           mask_phi(i)=0
2936           mask_theta(i)=0
2937          enddo
2938         endif
2939        enddo
2940        do j=1,nhfrag
2941         do i=hfrag(1,j),hfrag(2,j)
2942          mask(i)=0
2943          mask_phi(i)=0
2944          mask_theta(i)=0
2945         enddo
2946        enddo
2947        mask_r=.true.
2948
2949 !
2950 !      copy selected res from 1st to 2nd structure
2951 !
2952
2953        do i=1,nres
2954           if ( iff(i).eq.1 ) then
2955               theta(i)=theta1(i)
2956               phi(i)=phi1(i)
2957               alph(i)=alph1(i)
2958               omeg(i)=omeg1(i)
2959           endif
2960        enddo
2961
2962       if(debug) then
2963 !
2964 !     prepare description in linia variable
2965 !
2966         iwsk=0
2967         nf=0
2968         if (iff(1).eq.1) then
2969           iwsk=1
2970           nf=nf+1
2971           ij(nf)=1
2972         endif
2973         do i=2,nres
2974            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
2975              iwsk=1
2976              nf=nf+1
2977              ij(nf)=i
2978            endif
2979            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2980              iwsk=0
2981              nf=nf+1
2982              ij(nf)=i-1
2983            endif
2984         enddo
2985         if (iff(nres).eq.1) then
2986           nf=nf+1
2987           ij(nf)=nres
2988         endif
2989
2990         write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
2991                            "SELECT",ij(1)-1,"-",ij(2)-1,&
2992                                ",",ij(3)-1,"-",ij(4)-1
2993
2994       endif
2995 !
2996 !     run optimization
2997 !
2998       call contact_cp_min(var,ieval,in_pdb,linia,debug)
2999
3000       return
3001       end subroutine contact_cp
3002 !-----------------------------------------------------------------------------
3003       subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
3004
3005 !el      use minim
3006 !
3007 !    input : theta,phi,alph,omeg,in_pdb,linia,debug
3008 !    output : var,ieval
3009 !
3010 !      implicit real*8 (a-h,o-z)
3011 !      include 'DIMENSIONS'
3012       include 'mpif.h'
3013 !      include 'COMMON.SBRIDGE'
3014 !      include 'COMMON.FFIELD'
3015 !      include 'COMMON.IOUNITS'
3016 !      include 'COMMON.DISTFIT'
3017 !      include 'COMMON.VAR'
3018 !      include 'COMMON.CHAIN'
3019 !      include 'COMMON.MINIM'
3020
3021       character(len=50) :: linia
3022       integer :: nf,ij(4)
3023       real(kind=8) :: energy(0:n_ene)
3024       real(kind=8),dimension(6*nres) :: var     !(maxvar) (maxvar=6*maxres)
3025       real(kind=8) :: time0,time1
3026       integer :: ieval,info(3)
3027       logical :: debug,fail,reduce,change       !check_var,
3028
3029 !el local variables
3030       integer :: in_pdb,i,ipot0,ipot01,maxmin0,maxfun0,maxmin01,maxfun01,&
3031                  iretcode,nfun
3032       real(kind=8) :: wsc01,wscp01,welec01,wvdwpp01,wscloc01,wtor01,&
3033                  wtor_d01,wstrain0,etot
3034
3035        write(iout,'(a20,i6,a20)') &
3036                    '------------------',in_pdb,'-------------------'
3037 !el#ifdef MPI
3038        if (debug) then
3039         call chainbuild
3040         call write_pdb(1000+in_pdb,'combined structure',0d0)
3041 !el#ifdef MPI
3042         time0=MPI_WTIME()
3043 !el#endif
3044        endif
3045 !el#endif
3046 !
3047 !     run optimization of distances
3048 !
3049 !     uses d0(),w() and mask() for frozen 2D
3050 !
3051 !test---------------------------------------------
3052 !test       NX=NRES-3
3053 !test       NY=((NRES-4)*(NRES-5))/2
3054 !test       call distfit(debug,5000)
3055
3056        do i=1,nres
3057          mask_side(i)=0
3058        enddo
3059
3060        ipot01=ipot
3061        maxmin01=maxmin
3062        maxfun01=maxfun
3063 !       wstrain01=wstrain
3064        wsc01=wsc
3065        wscp01=wscp
3066        welec01=welec
3067        wvdwpp01=wvdwpp
3068 !      wang01=wang
3069        wscloc01=wscloc
3070        wtor01=wtor
3071        wtor_d01=wtor_d
3072
3073        ipot=6
3074        maxmin=2000
3075        maxfun=4000
3076 !       wstrain=1.0
3077        wsc=0.0
3078        wscp=0.0
3079        welec=0.0
3080        wvdwpp=0.0
3081 !      wang=0.0
3082        wscloc=0.0
3083        wtor=0.0
3084        wtor_d=0.0
3085
3086        call geom_to_var(nvar,var)
3087 !de       change=reduce(var)
3088        if (check_var(var,info)) then
3089           write(iout,*) 'cp_min error in input'
3090           print *,'cp_min error in input'
3091           return
3092        endif
3093
3094 !d       call etotal(energy(0))
3095 !d       call enerprint(energy(0))
3096 !d       call check_eint
3097 !el#ifdef MPI
3098        time0=MPI_WTIME()
3099 !dtest       call minimize(etot,var,iretcode,nfun)
3100 !dtest       write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
3101        time1=MPI_WTIME()
3102 !el#endif
3103 !d       call etotal(energy(0))
3104 !d       call enerprint(energy(0))
3105 !d       call check_eint
3106
3107        do i=1,nres
3108          mask_side(i)=1
3109        enddo
3110
3111        ipot=ipot01
3112        maxmin=maxmin01
3113        maxfun=maxfun01
3114 !       wstrain=wstrain01
3115        wsc=wsc01
3116        wscp=wscp01
3117        welec=welec01
3118        wvdwpp=wvdwpp01
3119 !      wang=wang01
3120        wscloc=wscloc01
3121        wtor=wtor01
3122        wtor_d=wtor_d01
3123 !test--------------------------------------------------
3124
3125        if(debug) then
3126 !el#ifdef MPI
3127         time1=MPI_WTIME()
3128 !el#endif
3129         write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
3130         call write_pdb(2000+in_pdb,'distfit structure',0d0)
3131        endif
3132
3133        ipot0=ipot
3134        maxmin0=maxmin
3135        maxfun0=maxfun
3136        wstrain0=wstrain
3137 !
3138 !      run soft pot. optimization
3139 !         with constrains:
3140 !             nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
3141 !         and frozen 2D:
3142 !             mask_phi(),mask_theta(),mask_side(),mask_r
3143 !
3144        ipot=6
3145        maxmin=2000
3146        maxfun=4000
3147 !el#ifdef MPI
3148 !de       change=reduce(var)
3149 !de       if (check_var(var,info)) write(iout,*) 'error before soft'
3150 !el#ifdef MPI
3151        time0=MPI_WTIME()
3152 !el#endif
3153        call minimize(etot,var,iretcode,nfun)
3154
3155        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3156 !el#ifdef MPI
3157        time1=MPI_WTIME()
3158 !el#endif
3159        write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,&
3160                nfun/(time1-time0),' SOFT eval/s'
3161
3162        if (debug) then
3163          call var_to_geom(nvar,var)
3164          call chainbuild
3165          call write_pdb(3000+in_pdb,'soft structure',etot)
3166        endif
3167 !el#endif
3168 !
3169 !      run full UNRES optimization with constrains and frozen 2D
3170 !      the same variables as soft pot. optimizatio
3171 !
3172        ipot=ipot0
3173        maxmin=maxmin0
3174        maxfun=maxfun0
3175 !
3176 ! check overlaps before calling full UNRES minim
3177 !
3178        call var_to_geom(nvar,var)
3179        call chainbuild
3180        call etotal(energy)
3181 #ifdef OSF
3182        write(iout,*) 'N7 ',energy(0)
3183        if (energy(0).ne.energy(0)) then
3184         write(iout,*) 'N7 error - gives NaN',energy(0)
3185        endif
3186 #endif
3187        ieval=1
3188        if (energy(1).eq.1.0d20) then
3189          write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
3190          call overlap_sc(fail)
3191          if(.not.fail) then
3192            call etotal(energy)
3193            ieval=ieval+1
3194            write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
3195          else
3196            mask_r=.false.
3197            nhpb= nhpb0
3198            link_start=1
3199            link_end=nhpb
3200            wstrain=wstrain0
3201            return
3202          endif
3203        endif
3204        call flush(iout)
3205 !
3206 !dte       time0=MPI_WTIME()
3207 !de       change=reduce(var)
3208 !de       if (check_var(var,info)) then
3209 !de         write(iout,*) 'error before mask dist'
3210 !de         call var_to_geom(nvar,var)
3211 !de         call chainbuild
3212 !de         call write_pdb(10000+in_pdb,'before mask dist',etot)
3213 !de       endif
3214 !dte       call minimize(etot,var,iretcode,nfun)
3215 !dte       write(iout,*)'SUMSL MASK DIST return code is',iretcode,
3216 !dte     &                          ' eval ',nfun
3217 !dte       ieval=ieval+nfun
3218 !dte
3219 !dte       time1=MPI_WTIME()
3220 !dte       write (iout,'(a,f6.2,f8.2,a)')
3221 !dte     &        ' Time for mask dist min.',time1-time0,
3222 !dte     &         nfun/(time1-time0),'  eval/s'
3223 !dte       call flush(iout)
3224 !el#ifdef MPI
3225        if (debug) then
3226          call var_to_geom(nvar,var)
3227          call chainbuild
3228          call write_pdb(4000+in_pdb,'mask dist',etot)
3229        endif
3230 !
3231 !      switch off freezing of 2D and
3232 !      run full UNRES optimization with constrains
3233 !
3234        mask_r=.false.
3235 !el#ifdef MPI
3236        time0=MPI_WTIME()
3237 !el#endif
3238 !de       change=reduce(var)
3239 !de       if (check_var(var,info)) then
3240 !de         write(iout,*) 'error before dist'
3241 !de         call var_to_geom(nvar,var)
3242 !de         call chainbuild
3243 !de         call write_pdb(11000+in_pdb,'before dist',etot)
3244 !de       endif
3245
3246        call minimize(etot,var,iretcode,nfun)
3247
3248 !de        change=reduce(var)
3249 !de        if (check_var(var,info)) then
3250 !de          write(iout,*) 'error after dist',ico
3251 !de          call var_to_geom(nvar,var)
3252 !de          call chainbuild
3253 !de          call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
3254 !de        endif
3255        write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
3256        ieval=ieval+nfun
3257
3258 !el#ifdef MPI
3259        time1=MPI_WTIME()
3260 !el#endif
3261        write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,&
3262                nfun/(time1-time0),'  eval/s'
3263
3264 !de       call etotal(energy(0))
3265 !de       write(iout,*) 'N7 after dist',energy(0)
3266        call flush(iout)
3267
3268        if (debug) then
3269         call var_to_geom(nvar,var)
3270         call chainbuild
3271         call write_pdb(in_pdb,linia,etot)
3272        endif
3273 !el#endif
3274 !
3275 !      reset constrains
3276 !
3277        nhpb= nhpb0
3278        link_start=1
3279        link_end=nhpb
3280        wstrain=wstrain0
3281
3282       return
3283       end subroutine contact_cp_min
3284 !-----------------------------------------------------------------------------
3285       subroutine softreg
3286
3287       use geometry, only:dist
3288 !el      use minim
3289 !      implicit real*8 (a-h,o-z)
3290 !      include 'DIMENSIONS'
3291       include 'mpif.h'
3292 !      include 'COMMON.GEO'
3293 !      include 'COMMON.CHAIN'
3294 !      include 'COMMON.IOUNITS'
3295 !      include 'COMMON.VAR'
3296 !      include 'COMMON.CONTROL'
3297 !      include 'COMMON.SBRIDGE'
3298 !      include 'COMMON.FFIELD'
3299 !      include 'COMMON.MINIM'
3300 !      include 'COMMON.INTERACT'
3301 !
3302 !      include 'COMMON.DISTFIT'
3303       integer :: iff(nres)
3304       real(kind=8) :: time0,time1
3305       real(kind=8) :: energy(0:n_ene),ee
3306       real(kind=8),dimension(6*nres) :: var     !(maxvar) (maxvar=6*maxres)
3307 !
3308       logical :: debug,ltest,fail
3309       character(len=50) :: linia
3310       integer :: ieval,i,j,in_pdb,ipot0,maxmin0,maxfun0,ico,nhpb_c,&
3311                  iretcode,nfun
3312       real(kind=8) :: wstrain0,wang0,etot
3313 !
3314       linia='test'
3315       debug=.true.
3316       in_pdb=0
3317
3318 !------------------------
3319 !
3320 !  freeze sec.elements
3321 !
3322        do i=1,nres
3323          mask_phi(i)=1
3324          mask_theta(i)=1
3325          mask_side(i)=1
3326          iff(i)=0
3327        enddo
3328
3329        do j=1,nbfrag
3330         do i=bfrag(1,j),bfrag(2,j)
3331          mask_phi(i)=0
3332          mask_theta(i)=0
3333          iff(i)=1
3334         enddo
3335         if (bfrag(3,j).le.bfrag(4,j)) then
3336          do i=bfrag(3,j),bfrag(4,j)
3337           mask_phi(i)=0
3338           mask_theta(i)=0
3339           iff(i)=1
3340          enddo
3341         else
3342          do i=bfrag(4,j),bfrag(3,j)
3343           mask_phi(i)=0
3344           mask_theta(i)=0
3345           iff(i)=1
3346          enddo
3347         endif
3348        enddo
3349        do j=1,nhfrag
3350         do i=hfrag(1,j),hfrag(2,j)
3351          mask_phi(i)=0
3352          mask_theta(i)=0
3353          iff(i)=1
3354         enddo
3355        enddo
3356        mask_r=.true.
3357
3358
3359
3360        nhpb0=nhpb
3361 !
3362 ! store dist. constrains
3363 !
3364        do i=1,nres-3
3365          do j=i+3,nres
3366            if ( iff(i).eq.1.and.iff(j).eq.1 ) then
3367             nhpb=nhpb+1
3368             ihpb(nhpb)=i
3369             jhpb(nhpb)=j
3370             forcon(nhpb)=0.1
3371             dhpb(nhpb)=DIST(i,j)
3372            endif
3373          enddo
3374        enddo
3375        call hpb_partition
3376
3377        if (debug) then
3378         call chainbuild
3379         call write_pdb(100+in_pdb,'input reg. structure',0d0)
3380        endif
3381
3382
3383        ipot0=ipot
3384        maxmin0=maxmin
3385        maxfun0=maxfun
3386        wstrain0=wstrain
3387        wang0=wang
3388 !
3389 !      run soft pot. optimization
3390 !
3391        ipot=6
3392        wang=3.0
3393        maxmin=2000
3394        maxfun=4000
3395        call geom_to_var(nvar,var)
3396 !el#ifdef MPI
3397        time0=MPI_WTIME()
3398 !el#endif
3399        call minimize(etot,var,iretcode,nfun)
3400
3401        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
3402 !el#ifdef MPI
3403        time1=MPI_WTIME()
3404 !el#endif
3405        write (iout,'(a,f6.2,f8.2,a)')'  Time for soft min.',time1-time0,&
3406                nfun/(time1-time0),' SOFT eval/s'
3407        if (debug) then
3408          call var_to_geom(nvar,var)
3409          call chainbuild
3410          call write_pdb(300+in_pdb,'soft structure',etot)
3411        endif
3412 !
3413 !      run full UNRES optimization with constrains and frozen 2D
3414 !      the same variables as soft pot. optimizatio
3415 !
3416        ipot=ipot0
3417        wang=wang0
3418        maxmin=maxmin0
3419        maxfun=maxfun0
3420 !el#ifdef MPI
3421        time0=MPI_WTIME()
3422 !el#endif
3423        call minimize(etot,var,iretcode,nfun)
3424        write(iout,*)'SUMSL MASK DIST return code is',iretcode,&
3425                                 ' eval ',nfun
3426        ieval=nfun
3427
3428 !el#ifdef MPI
3429        time1=MPI_WTIME()
3430 !el#endif
3431        write (iout,'(a,f6.2,f8.2,a)') &
3432               '  Time for mask dist min.',time1-time0,&
3433                nfun/(time1-time0),'  eval/s'
3434        if (debug) then
3435          call var_to_geom(nvar,var)
3436          call chainbuild
3437          call write_pdb(400+in_pdb,'mask & dist',etot)
3438        endif
3439 !
3440 !      switch off constrains and
3441 !      run full UNRES optimization with frozen 2D
3442 !
3443
3444 !
3445 !      reset constrains
3446 !
3447        nhpb_c=nhpb
3448        nhpb=nhpb0
3449        link_start=1
3450        link_end=nhpb
3451        wstrain=wstrain0
3452
3453 !el#ifdef MPI
3454        time0=MPI_WTIME()
3455 !el#endif
3456        call minimize(etot,var,iretcode,nfun)
3457        write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
3458        ieval=ieval+nfun
3459
3460 !el#ifdef MPI
3461        time1=MPI_WTIME()
3462 !el#endif
3463        write (iout,'(a,f6.2,f8.2,a)')'  Time for mask min.',time1-time0,&
3464                nfun/(time1-time0),'  eval/s'
3465
3466
3467        if (debug) then
3468         call var_to_geom(nvar,var)
3469         call chainbuild
3470         call write_pdb(500+in_pdb,'mask 2d frozen',etot)
3471        endif
3472
3473        mask_r=.false.
3474
3475
3476 !
3477 !      run full UNRES optimization with constrains and NO frozen 2D
3478 !
3479
3480        nhpb=nhpb_c
3481        link_start=1
3482        link_end=nhpb
3483        maxfun=maxfun0/5
3484
3485        do ico=1,5
3486
3487        wstrain=wstrain0/ico
3488
3489 !el#ifdef MPI
3490        time0=MPI_WTIME()
3491 !el#endif
3492        call minimize(etot,var,iretcode,nfun)
3493        write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3494          ' SUMSL DIST',wstrain,' return code is',iretcode,&
3495                                 ' eval ',nfun
3496        ieval=nfun
3497
3498 !el#ifdef MPI
3499        time1=MPI_WTIME()
3500 !el#endif
3501        write (iout,'(a,f6.2,f8.2,a)') &
3502               '  Time for dist min.',time1-time0,&
3503                nfun/(time1-time0),'  eval/s'
3504        if (debug) then
3505          call var_to_geom(nvar,var)
3506          call chainbuild
3507          call write_pdb(600+in_pdb+ico,'dist cons',etot)
3508        endif
3509
3510        enddo
3511 !
3512        nhpb=nhpb0
3513        link_start=1
3514        link_end=nhpb
3515        wstrain=wstrain0
3516        maxfun=maxfun0
3517
3518
3519 !
3520       if (minim) then
3521 !el#ifdef MPI
3522        time0=MPI_WTIME()
3523 !el#endif
3524        call minimize(etot,var,iretcode,nfun)
3525        write(iout,*)'------------------------------------------------'
3526        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,&
3527         '+ DIST eval',ieval
3528
3529 !el#ifdef MPI
3530        time1=MPI_WTIME()
3531 !el#endif
3532        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,&
3533                nfun/(time1-time0),' eval/s'
3534
3535
3536        call var_to_geom(nvar,var)
3537        call chainbuild
3538        call write_pdb(999,'full min',etot)
3539       endif
3540 !el#endif
3541       return
3542       end subroutine softreg
3543 !-----------------------------------------------------------------------------
3544       subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
3545
3546       use geometry, only:dist
3547 !el      use minim
3548 !      implicit real*8 (a-h,o-z)
3549 !      include 'DIMENSIONS'
3550       include 'mpif.h'
3551 !      include 'COMMON.GEO'
3552 !      include 'COMMON.VAR'
3553 !      include 'COMMON.INTERACT'
3554 !      include 'COMMON.IOUNITS'
3555 !      include 'COMMON.DISTFIT'
3556 !      include 'COMMON.SBRIDGE'
3557 !      include 'COMMON.CONTROL'
3558 !      include 'COMMON.FFIELD'
3559 !      include 'COMMON.MINIM'
3560 !      include 'COMMON.CHAIN'
3561       real(kind=8) :: time0,time1
3562       real(kind=8) :: energy(0:n_ene),ee
3563       real(kind=8),dimension(6*nres) :: var     !(maxvar) (maxvar=6*maxres)
3564       integer :: jdata(5),isec(nres)
3565 !
3566 !el local variables
3567       integer :: i1,i2,i3,i4,i5,ieval,ij
3568       integer :: i,j,nft_sc,ishift,iretcode,nfun,maxfun0,ico
3569       real(kind=8) :: etot,wscloc0,wstrain0
3570
3571       jdata(1)=i1
3572       jdata(2)=i2
3573       jdata(3)=i3
3574       jdata(4)=i4
3575       jdata(5)=i5
3576
3577       call secondary2(.false.)
3578
3579       do i=1,nres
3580           isec(i)=0
3581       enddo
3582       do j=1,nbfrag
3583        do i=bfrag(1,j),bfrag(2,j)
3584           isec(i)=1
3585        enddo
3586        do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
3587           isec(i)=1
3588        enddo
3589       enddo
3590       do j=1,nhfrag
3591        do i=hfrag(1,j),hfrag(2,j)
3592           isec(i)=2
3593        enddo
3594       enddo
3595
3596 !
3597 ! cut strands at the ends
3598 !
3599       if (jdata(2)-jdata(1).gt.3) then
3600        jdata(1)=jdata(1)+1
3601        jdata(2)=jdata(2)-1
3602        if (jdata(3).lt.jdata(4)) then
3603           jdata(3)=jdata(3)+1
3604           jdata(4)=jdata(4)-1
3605        else
3606           jdata(3)=jdata(3)-1
3607           jdata(4)=jdata(4)+1
3608        endif
3609       endif
3610
3611 !v      call chainbuild
3612 !v      call etotal(energy(0))
3613 !v      etot=energy(0)
3614 !v      write(iout,*) nnt,nct,etot
3615 !v      call write_pdb(ij*100,'first structure',etot)
3616 !v      write(iout,*) 'N16 test',(jdata(i),i=1,5)
3617
3618 !------------------------
3619 !      generate constrains
3620 !
3621        ishift=jdata(5)-2
3622        if(ishift.eq.0) ishift=-2
3623        nhpb0=nhpb
3624        call chainbuild
3625        do i=jdata(1),jdata(2)
3626         isec(i)=-1
3627         if(jdata(4).gt.jdata(3))then
3628          do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
3629             isec(j)=-1
3630 !d            print *,i,j,j+ishift
3631             nhpb=nhpb+1
3632             ihpb(nhpb)=i
3633             jhpb(nhpb)=j
3634             forcon(nhpb)=1000.0
3635             dhpb(nhpb)=DIST(i,j+ishift)
3636          enddo
3637         else
3638          do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
3639             isec(j)=-1
3640 !d            print *,i,j,j+ishift
3641             nhpb=nhpb+1
3642             ihpb(nhpb)=i
3643             jhpb(nhpb)=j
3644             forcon(nhpb)=1000.0
3645             dhpb(nhpb)=DIST(i,j+ishift)
3646          enddo
3647         endif
3648        enddo
3649
3650        do i=nnt,nct-2
3651          do j=i+2,nct
3652            if(isec(i).gt.0.or.isec(j).gt.0) then
3653 !d            print *,i,j
3654             nhpb=nhpb+1
3655             ihpb(nhpb)=i
3656             jhpb(nhpb)=j
3657             forcon(nhpb)=0.1
3658             dhpb(nhpb)=DIST(i,j)
3659            endif
3660          enddo
3661        enddo
3662
3663        call hpb_partition
3664
3665        call geom_to_var(nvar,var)
3666        maxfun0=maxfun
3667        wstrain0=wstrain
3668        maxfun=4000/5
3669
3670        do ico=1,5
3671
3672         wstrain=wstrain0/ico
3673
3674 !v        time0=MPI_WTIME()
3675         call minimize(etot,var,iretcode,nfun)
3676         write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3677          ' SUMSL DIST',wstrain,' return code is',iretcode,&
3678                                 ' eval ',nfun
3679         ieval=ieval+nfun
3680 !v        time1=MPI_WTIME()
3681 !v       write (iout,'(a,f6.2,f8.2,a)')
3682 !v     &        '  Time for dist min.',time1-time0,
3683 !v     &         nfun/(time1-time0),'  eval/s'
3684 !v         call var_to_geom(nvar,var)
3685 !v         call chainbuild
3686 !v         call write_pdb(ij*100+ico,'dist cons',etot)
3687
3688        enddo
3689 !
3690        nhpb=nhpb0
3691        call hpb_partition
3692        wstrain=wstrain0
3693        maxfun=maxfun0
3694 !
3695 !d      print *,etot
3696       wscloc0=wscloc
3697       wscloc=10.0
3698       call sc_move(nnt,nct,100,100d0,nft_sc,etot)
3699       wscloc=wscloc0
3700 !v      call chainbuild
3701 !v      call etotal(energy(0))
3702 !v      etot=energy(0)
3703 !v      call write_pdb(ij*100+10,'sc_move',etot)
3704 !d      call intout
3705 !d      print *,nft_sc,etot
3706
3707       return
3708       end subroutine beta_slide
3709 !-----------------------------------------------------------------------------
3710       subroutine beta_zip(i1,i2,ieval,ij)
3711
3712 !el      use minim
3713 !      implicit real*8 (a-h,o-z)
3714 !      include 'DIMENSIONS'
3715       include 'mpif.h'
3716 !      include 'COMMON.GEO'
3717 !      include 'COMMON.VAR'
3718 !      include 'COMMON.INTERACT'
3719 !      include 'COMMON.IOUNITS'
3720 !      include 'COMMON.DISTFIT'
3721 !      include 'COMMON.SBRIDGE'
3722 !      include 'COMMON.CONTROL'
3723 !      include 'COMMON.FFIELD'
3724 !      include 'COMMON.MINIM'
3725 !      include 'COMMON.CHAIN'
3726       real(kind=8) :: time0,time1
3727       real(kind=8) :: energy(0:n_ene),ee
3728       real(kind=8),dimension(6*nres) :: var     !(maxvar) (maxvar=6*maxres)
3729       character(len=10) :: test
3730 !el local variables
3731       integer :: i1,i2,ieval,ij,ico,iretcode,nfun,maxfun0
3732       real(kind=8) :: etot,wstrain0
3733 !v      call chainbuild
3734 !v      call etotal(energy(0))
3735 !v      etot=energy(0)
3736 !v      write(test,'(2i5)') i1,i2
3737 !v      call write_pdb(ij*100,test,etot)
3738 !v      write(iout,*) 'N17 test',i1,i2,etot,ij
3739
3740 !
3741 !      generate constrains
3742 !
3743        nhpb0=nhpb
3744        nhpb=nhpb+1
3745        ihpb(nhpb)=i1
3746        jhpb(nhpb)=i2
3747        forcon(nhpb)=1000.0
3748        dhpb(nhpb)=4.0
3749
3750        call hpb_partition
3751
3752        call geom_to_var(nvar,var)
3753        maxfun0=maxfun
3754        wstrain0=wstrain
3755        maxfun=1000/5
3756
3757        do ico=1,5
3758         wstrain=wstrain0/ico
3759 !v        time0=MPI_WTIME()
3760         call minimize(etot,var,iretcode,nfun)
3761         write(iout,'(a10,f6.3,a14,i3,a6,i5)') &
3762          ' SUMSL DIST',wstrain,' return code is',iretcode,&
3763                                 ' eval ',nfun
3764         ieval=ieval+nfun
3765 !v        time1=MPI_WTIME()
3766 !v       write (iout,'(a,f6.2,f8.2,a)')
3767 !v     &        '  Time for dist min.',time1-time0,
3768 !v     &         nfun/(time1-time0),'  eval/s'
3769 ! do not comment the next line
3770          call var_to_geom(nvar,var)
3771 !v         call chainbuild
3772 !v         call write_pdb(ij*100+ico,'dist cons',etot)
3773        enddo
3774
3775        nhpb=nhpb0
3776        call hpb_partition
3777        wstrain=wstrain0
3778        maxfun=maxfun0
3779
3780 !v      call etotal(energy(0))
3781 !v      etot=energy(0)
3782 !v      write(iout,*) 'N17 test end',i1,i2,etot,ij
3783
3784       return
3785       end subroutine beta_zip
3786 !-----------------------------------------------------------------------------
3787 ! thread.F
3788 !-----------------------------------------------------------------------------
3789       subroutine thread_seq
3790
3791       use geometry, only:dist
3792       use random, only:iran_num
3793       use control, only:tcpu
3794       use regularize_, only:regularize
3795       use mcm_data, only: nsave_part,nacc_tot
3796 ! Thread the sequence through a database of known structures
3797 !      implicit real*8 (a-h,o-z)
3798 !      include 'DIMENSIONS'
3799       use MPI_data      !include 'COMMON.INFO'
3800       use MPI_
3801 #ifdef MPI
3802       include 'mpif.h'
3803 #endif
3804 !      include 'COMMON.CONTROL'
3805 !      include 'COMMON.CHAIN'
3806 !      include 'COMMON.DBASE'
3807 !      include 'COMMON.INTERACT'
3808 !      include 'COMMON.VAR'
3809 !      include 'COMMON.THREAD'
3810 !      include 'COMMON.FFIELD'
3811 !      include 'COMMON.SBRIDGE'
3812 !      include 'COMMON.HEADER'
3813 !      include 'COMMON.IOUNITS'
3814 !      include 'COMMON.TIME1'
3815 !      include 'COMMON.CONTACTS'
3816 !      include 'COMMON.MCM'
3817 !      include 'COMMON.NAMES'
3818 #ifdef MPI
3819       integer :: ThreadId,ThreadType,Kwita
3820 #endif
3821       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
3822       real(kind=8) :: przes(3),obr(3,3)
3823       real(kind=8) :: time_for_thread
3824       logical :: found_pattern,non_conv
3825       character(len=32) :: head_pdb
3826       real(kind=8) :: energia(0:n_ene)
3827       integer :: i,j,ithread,itrial,ii,jj,nres_t,ist,ipattern,iretcode,&
3828             link_end0,iproc
3829       real(kind=8) :: dcj,rms,frac,frac_nn,co,etot,curr_tim,curr_tim1
3830
3831       n_ene_comp=nprint_ene
3832 !
3833 ! Body
3834 !
3835 #ifdef MPI
3836       if (me.eq.king) then
3837         do i=1,nctasks
3838           nsave_part(i)=0
3839         enddo
3840       endif
3841       nacc_tot=0
3842
3843       Kwita=0
3844 #endif
3845       close(igeom)
3846       close(ipdb)
3847       close(istat)
3848       do i=1,maxthread
3849         do j=1,14
3850           ener0(j,i)=0.0D0
3851           ener(j,i)=0.0D0
3852         enddo
3853       enddo
3854       nres0=nct-nnt+1
3855       ave_time_for_thread=0.0D0
3856       max_time_for_thread=0.0D0
3857 !d    print *,'nthread=',nthread,' nseq=',nseq,' nres0=',nres0
3858       nthread=nexcl+nthread
3859       do ithread=1,nthread
3860         found_pattern=.false.
3861         itrial=0
3862         do while (.not.found_pattern)
3863           itrial=itrial+1
3864           if (itrial.gt.1000) then
3865             write (iout,'(/a/)') 'Too many attempts to find pattern.'
3866             nthread=ithread-1
3867 #ifdef MPI
3868             call recv_stop_sig(Kwita)
3869             call send_stop_sig(-3)
3870 #endif
3871             goto 777
3872           endif
3873 ! Find long enough chain in the database
3874           ii=iran_num(1,nseq)
3875           nres_t=nres_base(1,ii)
3876 ! Select the starting position to thread.
3877           print *,'nseq',nseq,' ii=',ii,' nres_t=',&
3878             nres_t,' nres0=',nres0
3879           if (nres_t.ge.nres0) then
3880             ist=iran_num(0,nres_t-nres0)
3881 #ifdef MPI
3882             if (Kwita.eq.0) call recv_stop_sig(Kwita)
3883             if (Kwita.lt.0) then
3884               write (iout,*) 'Stop signal received. Terminating.'
3885               write (*,*) 'Stop signal received. Terminating.'
3886               nthread=ithread-1
3887               write (*,*) 'ithread=',ithread,' nthread=',nthread
3888               goto 777
3889             endif
3890             call pattern_receive
3891 #endif
3892             do i=1,nexcl
3893               if (iexam(1,i).eq.ii .and. iexam(2,i).eq.ist) goto 10
3894             enddo
3895             found_pattern=.true.
3896           endif
3897 ! If this point is reached, the pattern has not yet been examined.
3898    10     continue
3899 !         print *,'found_pattern:',found_pattern
3900         enddo
3901         nexcl=nexcl+1
3902         iexam(1,nexcl)=ii
3903         iexam(2,nexcl)=ist
3904 #ifdef MPI
3905         if (Kwita.eq.0) call recv_stop_sig(Kwita)
3906         if (Kwita.lt.0) then
3907           write (iout,*) 'Stop signal received. Terminating.'
3908           nthread=ithread-1
3909           write (*,*) 'ithread=',ithread,' nthread=',nthread
3910           goto 777
3911         endif
3912         call pattern_send
3913 #endif
3914         ipatt(1,ithread)=ii
3915         ipatt(2,ithread)=ist
3916 #ifdef MPI
3917         write (iout,'(/80(1h*)/a,i4,a,i5,2a,i3,a,i3,a,i3/)') &
3918          'Processor:',me,' Attempt:',ithread,&
3919          ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3920          ' start at res.',ist+1
3921         write (*,'(a,i4,a,i5,2a,i3,a,i3,a,i3)') 'Processor:',me,&
3922          ' Attempt:',ithread,&
3923          ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3924          ' start at res.',ist+1
3925 #else
3926         write (iout,'(/80(1h*)/a,i5,2a,i3,a,i3,a,i3/)') &
3927          'Attempt:',ithread,&
3928          ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3929          ' start at res.',ist+1
3930         write (*,'(a,i5,2a,i3,a,i3,a,i3)') &
3931          'Attempt:',ithread,&
3932          ' pattern: ',str_nam(ii),nres_base(2,ii),':',nres_base(3,ii),&
3933          ' start at res.',ist+1
3934 #endif
3935         ipattern=ii
3936 ! Copy coordinates from the database.
3937         ist=ist-(nnt-1)
3938         do i=nnt,nct
3939           do j=1,3
3940             c(j,i)=cart_base(j,i+ist,ii)
3941 !           cref(j,i)=c(j,i)
3942           enddo
3943 !d        write (iout,'(a,i4,3f10.5)') restyp(itype(i,1)),i,(c(j,i),j=1,3)
3944         enddo
3945 !d      call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,
3946 !d             non_conv)
3947 !d      write (iout,'(a,f10.5)')
3948 !d   &  'Initial RMS deviation from reference structure:',rms
3949         if (itype(nres,1).eq.ntyp1) then
3950           do j=1,3
3951             dcj=c(j,nres-2)-c(j,nres-3)
3952             c(j,nres)=c(j,nres-1)+dcj
3953             c(j,2*nres)=c(j,nres)
3954           enddo
3955         endif
3956         if (itype(1,1).eq.ntyp1) then
3957           do j=1,3
3958             dcj=c(j,4)-c(j,3)
3959             c(j,1)=c(j,2)-dcj
3960             c(j,nres+1)=c(j,1)
3961           enddo
3962         endif
3963         call int_from_cart(.false.,.false.)
3964 !d      print *,'Exit INT_FROM_CART.'
3965 !d      print *,'nhpb=',nhpb
3966         do i=nss+1,nhpb
3967           ii=ihpb(i)
3968           jj=jhpb(i)
3969           dhpb(i)=dist(ii,jj)
3970 !         write (iout,'(2i5,2f10.5)') ihpb(i),jhpb(i),dhpb(i),forcon(i)
3971         enddo
3972 !       stop 'End generate'
3973 ! Generate SC conformations.
3974         call sc_conf
3975 !       call intout
3976 #ifdef MPI
3977 !d      print *,'Processor:',me,': exit GEN_SIDE.'
3978 #else
3979 !d      print *,'Exit GEN_SIDE.'
3980 #endif
3981 ! Calculate initial energy.
3982         call chainbuild
3983         call etotal(energia)
3984         etot=energia(0)
3985         do i=1,n_ene_comp
3986           ener0(i,ithread)=energia(i)
3987         enddo
3988         ener0(n_ene_comp+1,ithread)=energia(0)
3989         if (refstr) then
3990           call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
3991           ener0(n_ene_comp+3,ithread)=contact_fract(ncont,ncont_ref,&
3992               icont,icont_ref)
3993           ener0(n_ene_comp+2,ithread)=rms
3994           ener0(n_ene_comp+4,ithread)=frac
3995           ener0(n_ene_comp+5,ithread)=frac_nn
3996         endif
3997         ener0(n_ene_comp+3,ithread)=0.0d0
3998 ! Minimize energy.
3999 #ifdef MPI
4000        print*,'Processor:',me,' ithread=',ithread,' Start REGULARIZE.'
4001 #else
4002         print*,'ithread=',ithread,' Start REGULARIZE.'
4003 #endif
4004         curr_tim=tcpu()
4005         call regularize(nct-nnt+1,etot,rms,&
4006                         cart_base(1,ist+nnt,ipattern),iretcode)
4007         curr_tim1=tcpu()
4008         time_for_thread=curr_tim1-curr_tim
4009         ave_time_for_thread= &
4010         ((ithread-1)*ave_time_for_thread+time_for_thread)/ithread
4011         if (time_for_thread.gt.max_time_for_thread) &
4012          max_time_for_thread=time_for_thread
4013 #ifdef MPI
4014         print *,'Processor',me,': Exit REGULARIZE.'
4015         if (WhatsUp.eq.2) then
4016           write (iout,*) &
4017         'Sufficient number of confs. collected. Terminating.'
4018           nthread=ithread-1
4019           goto 777
4020         else if (WhatsUp.eq.-1) then
4021           nthread=ithread-1
4022           write (iout,*) 'Time up in REGULARIZE. Call SEND_STOP_SIG.'
4023           if (Kwita.eq.0) call recv_stop_sig(Kwita)
4024           call send_stop_sig(-2)
4025           goto 777
4026         else if (WhatsUp.eq.-2) then
4027           nthread=ithread-1
4028           write (iout,*) 'Timeup signal received. Terminating.'
4029           goto 777
4030         else if (WhatsUp.eq.-3) then
4031           nthread=ithread-1
4032           write (iout,*) 'Error stop signal received. Terminating.'
4033           goto 777
4034         endif
4035 #else
4036         print *,'Exit REGULARIZE.'
4037         if (iretcode.eq.11) then
4038           write (iout,'(/a/)') &
4039       '******* Allocated time exceeded in SUMSL. The program will stop.'
4040           nthread=ithread-1
4041           goto 777
4042         endif
4043 #endif
4044         head_pdb=titel(:24)//':'//str_nam(ipattern)
4045         if (outpdb) call pdbout(etot,head_pdb,ipdb)
4046         if (outmol2) call mol2out(etot,head_pdb)
4047 !       call intout
4048         call briefout(ithread,etot)
4049         link_end0=link_end
4050         link_end=min0(link_end,nss)
4051         write (iout,*) 'link_end=',link_end,' link_end0=',link_end0,&
4052                        ' nss=',nss
4053         call etotal(energia)
4054 !       call enerprint(energia(0))
4055         link_end=link_end0
4056 !d      call chainbuild
4057 !d      call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obr,non_conv)
4058 !d      write (iout,'(a,f10.5)')
4059 !d   &  'RMS deviation from reference structure:',dsqrt(rms)
4060         do i=1,n_ene_comp
4061           ener(i,ithread)=energia(i)
4062         enddo
4063         ener(n_ene_comp+1,ithread)=energia(0)
4064         ener(n_ene_comp+3,ithread)=rms
4065         if (refstr) then
4066           call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
4067           ener(n_ene_comp+2,ithread)=rms
4068           ener(n_ene_comp+4,ithread)=frac
4069           ener(n_ene_comp+5,ithread)=frac_nn
4070         endif
4071         call write_stat_thread(ithread,ipattern,ist)
4072 !        write (istat,'(i4,2x,a8,i4,11(1pe14.5),2(0pf8.3),f8.5)')
4073 !     &  ithread,str_nam(ipattern),ist+1,(ener(k,ithread),k=1,11),
4074 !     &  (ener(k,ithread),k=12,14)
4075 #ifdef MPI
4076         if (me.eq.king) then
4077           nacc_tot=nacc_tot+1
4078           call pattern_receive
4079           call receive_MCM_info
4080           if (nacc_tot.ge.nthread) then
4081             write (iout,*) &
4082            'Sufficient number of conformations collected nacc_tot=',&
4083            nacc_tot,'. Stopping other processors and terminating.'
4084             write (*,*) &
4085            'Sufficient number of conformations collected nacc_tot=',&
4086            nacc_tot,'. Stopping other processors and terminating.'
4087            call recv_stop_sig(Kwita)
4088            if (Kwita.eq.0) call send_stop_sig(-1)
4089            nthread=ithread
4090            goto 777
4091           endif
4092         else
4093           call send_MCM_info(2)
4094         endif
4095 #endif
4096         if (timlim-curr_tim1-safety .lt. max_time_for_thread) then
4097           write (iout,'(/2a)') &
4098        '********** There would be not enough time for another thread. ',&
4099        'The program will stop.'
4100           write (*,'(/2a)') &
4101        '********** There would be not enough time for another thread. ',&
4102        'The program will stop.'
4103           write (iout,'(a,1pe14.4/)') &
4104           'Elapsed time for last threading step: ',time_for_thread
4105           nthread=ithread
4106 #ifdef MPI
4107           call recv_stop_sig(Kwita)
4108           call send_stop_sig(-2)
4109 #endif
4110           goto 777
4111         else
4112           curr_tim=curr_tim1
4113           write (iout,'(a,1pe14.4)') &
4114           'Elapsed time for this threading step: ',time_for_thread
4115         endif
4116 #ifdef MPI
4117         if (Kwita.eq.0) call recv_stop_sig(Kwita)
4118         if (Kwita.lt.0) then
4119           write (iout,*) 'Stop signal received. Terminating.'
4120           write (*,*) 'Stop signal received. Terminating.'
4121           nthread=ithread
4122           write (*,*) 'nthread=',nthread,' ithread=',ithread
4123           goto 777
4124         endif
4125 #endif
4126       enddo
4127 #ifdef MPI
4128       call send_stop_sig(-1)
4129 #endif
4130   777 continue
4131 #ifdef MPI
4132 ! Any messages left for me?
4133       call pattern_receive
4134       if (Kwita.eq.0) call recv_stop_sig(Kwita)
4135 #endif
4136       call write_thread_summary
4137 #ifdef MPI
4138       if (king.eq.king) then
4139         Kwita=1
4140         do while (Kwita.ne.0 .or. nacc_tot.ne.0)
4141           Kwita=0
4142           nacc_tot=0
4143           call recv_stop_sig(Kwita)
4144           call receive_MCM_info
4145         enddo
4146         do iproc=1,nprocs-1
4147           call receive_thread_results(iproc)
4148         enddo
4149         call write_thread_summary
4150       else
4151         call send_thread_results
4152       endif
4153 #endif
4154       return
4155       end subroutine thread_seq
4156 !-----------------------------------------------------------------------------
4157       subroutine sc_conf
4158
4159 ! Sample (hopefully) optimal SC orientations given backcone conformation.
4160 !el      use comm_srutu
4161       use random, only:iran_num
4162 !      implicit real*8 (a-h,o-z)
4163 !      include 'DIMENSIONS'
4164 !      include 'COMMON.CHAIN'
4165 !      include 'COMMON.DBASE'
4166 !      include 'COMMON.INTERACT'
4167 !      include 'COMMON.VAR'
4168 !      include 'COMMON.THREAD'
4169 !      include 'COMMON.FFIELD'
4170 !      include 'COMMON.SBRIDGE'
4171 !      include 'COMMON.HEADER'
4172 !      include 'COMMON.GEO'
4173 !      include 'COMMON.IOUNITS'
4174       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
4175 !el      integer :: icall
4176 !el      common /srutu/ icall
4177       real(kind=8) :: energia(0:n_ene)
4178       logical :: glycine,fail
4179       integer :: i,maxsample,link_end0,ind_sc,isample
4180       real(kind=8) :: alph0,omeg0,e1,e0
4181
4182       maxsample=10
4183       link_end0=link_end
4184       link_end=min0(link_end,nss)
4185       do i=nnt,nct
4186         if (itype(i,1).ne.10) then
4187 !d        print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
4188           call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail)
4189         endif
4190       enddo
4191       call chainbuild
4192       call etotal(energia)
4193       e0 = energia(0)
4194       do isample=1,maxsample
4195 ! Choose a non-glycine side chain.
4196         glycine=.true.
4197         do while(glycine)
4198           ind_sc=iran_num(nnt,nct)
4199           glycine=(itype(ind_sc,1).eq.10)
4200         enddo
4201         alph0=alph(ind_sc)
4202         omeg0=omeg(ind_sc)
4203         call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),&
4204              omeg(ind_sc),fail)
4205         call chainbuild
4206         call etotal(energia)
4207 !d      write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')
4208 !d   &   'Step:',isample,' SC',ind_sc,' alpha',alph(ind_sc)*rad2deg,
4209 !d   &   ' omega',omeg(ind_sc)*rad2deg,' old energy',e0,' new energy',e1
4210         e1=energia(0)
4211         if (e0.le.e1) then
4212           alph(ind_sc)=alph0
4213           omeg(ind_sc)=omeg0
4214         else
4215           e0=e1
4216         endif
4217       enddo
4218       link_end=link_end0
4219       return
4220       end subroutine sc_conf
4221 !-----------------------------------------------------------------------------
4222 ! minim_jlee.F
4223 !-----------------------------------------------------------------------------
4224       logical function check_var(var,info)
4225
4226       use MPI_data
4227       use geometry_data
4228 !      implicit real*8 (a-h,o-z)
4229 !      include 'DIMENSIONS'
4230 !      include 'COMMON.VAR'
4231 !      include 'COMMON.IOUNITS'
4232 !      include 'COMMON.GEO'
4233 !      include 'COMMON.SETUP'
4234       real(kind=8),dimension(6*nres) :: var     !(maxvar) (maxvar=6*maxres)
4235       integer,dimension(3) :: info
4236       integer :: i,j
4237 ! AL -------
4238        check_var=.false.
4239        do i=nphi+ntheta+1,nphi+ntheta+nside
4240 ! Check the side chain "valence" angles alpha
4241          if (var(i).lt.1.0d-7) then
4242            write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4243            write (iout,*) 'Processor',me,'received bad variables!!!!'
4244            write (iout,*) 'Variables'
4245            write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4246            write (iout,*) 'Continuing calculations at this point',&
4247        ' could destroy the results obtained so far... ABORTING!!!!!!'
4248            write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4249            'valence angle alpha',i-nphi-ntheta,var(i),&
4250            'n it',info(1),info(2),'mv ',info(3)
4251            write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4252            write (*,*) 'Processor',me,'received bad variables!!!!'
4253            write (*,*) 'Variables'
4254            write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4255            write (*,*) 'Continuing calculations at this point',&
4256        ' could destroy the results obtained so far... ABORTING!!!!!!'
4257            write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4258            'valence angle alpha',i-nphi-ntheta,var(i),&
4259            'n it',info(1),info(2),'mv ',info(3)
4260            check_var=.true.
4261            return
4262          endif
4263        enddo
4264 ! Check the backbone "valence" angles theta
4265        do i=nphi+1,nphi+ntheta
4266          if (var(i).lt.1.0d-7) then
4267            write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4268            write (iout,*) 'Processor',me,'received bad variables!!!!'
4269            write (iout,*) 'Variables'
4270            write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4271            write (iout,*) 'Continuing calculations at this point',&
4272        ' could destroy the results obtained so far... ABORTING!!!!!!'
4273            write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4274            'valence angle theta',i-nphi,var(i),&
4275            'n it',info(1),info(2),'mv ',info(3)
4276            write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
4277            write (*,*) 'Processor',me,'received bad variables!!!!'
4278            write (*,*) 'Variables'
4279            write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
4280            write (*,*) 'Continuing calculations at this point',&
4281        ' could destroy the results obtained so far... ABORTING!!!!!!'
4282            write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') &
4283            'valence angle theta',i-nphi,var(i),&
4284            'n it',info(1),info(2),'mv ',info(3)
4285            check_var=.true.
4286            return
4287          endif
4288        enddo
4289       return
4290       end function check_var
4291 !-----------------------------------------------------------------------------
4292 ! distfit.f
4293 !-----------------------------------------------------------------------------
4294       subroutine distfit(debug,maxit)
4295
4296       use geometry_data, only: phi
4297       use compare_data
4298       use md_calc
4299 !      implicit real*8 (a-h,o-z)
4300 !      include 'DIMENSIONS'
4301 !      include 'COMMON.CHAIN'
4302 !      include 'COMMON.VAR'
4303 !      include 'COMMON.IOUNITS'
4304 !      include 'COMMON.DISTFIT'
4305       integer :: i,maxit,MAXMAR,IT,IMAR
4306       real(kind=8),DIMENSION(nres) :: X,DIAGH,phiold    !(maxres)
4307       logical :: debug,sing
4308       real(kind=8) :: TOL,RL,F0,AIN,F1
4309
4310 !input------------------------------------
4311 !       NX=NRES-3
4312 !       NY=((NRES-4)*(NRES-5))/2
4313 !input------------------------------------
4314 !test      MAXIT=20
4315       TOL=0.5
4316       MAXMAR=10
4317       RL=100.0
4318
4319       CALL TRANSFER(NRES,phi,phiold)
4320
4321       F0=RDIF()
4322
4323 !d      WRITE (IOUT,*) 'DISTFIT: F0=',F0
4324
4325
4326       DO IT=1,MAXIT
4327         CALL RDERIV
4328         CALL HEVAL
4329
4330         DO I=1,NX
4331           DIAGH(I)=H(I,I)
4332         ENDDO
4333         RL=RL*0.1
4334
4335         DO IMAR=1,MAXMAR
4336           DO I=1,NX
4337             H(I,I)=DIAGH(I)+RL
4338           ENDDO
4339           CALL TRANSFER(NX,XX,X)
4340           CALL BANACH(NX,NRES,H,X,sing)
4341           AIN=0.0
4342           DO I=1,NX
4343             AIN=AIN+DABS(X(I))
4344           ENDDO
4345           IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN
4346               if (debug) then
4347               WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED'
4348               WRITE (IOUT,*) 'IT=',it,'F=',F0
4349               endif
4350             RETURN
4351           ENDIF
4352           DO I=4,NRES
4353             phi(I)=phiold(I)+mask(i)*X(I-3)
4354 !            print *,X(I-3)
4355           ENDDO
4356
4357           F1=RDIF()
4358 !d          WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1
4359           IF (F1.LT.F0) THEN
4360             CALL TRANSFER(NRES,phi,phiold)
4361             F0=F1
4362             GOTO 1
4363           ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN
4364             if (debug) then
4365             WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT'
4366             WRITE (IOUT,*) 'IT=',it,'F=',F1
4367             endif
4368             RETURN
4369           ENDIF
4370           RL=RL*10.0
4371         ENDDO
4372         WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED'
4373         WRITE (IOUT,*) 'IT=',it,'F=',F0
4374         CALL TRANSFER(NRES,phiold,phi)
4375         RETURN
4376     1   continue
4377 !d        write (iout,*) "it",it," imar",imar," f0",f0
4378       enddo
4379       WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit
4380       return
4381       end subroutine distfit
4382 !-----------------------------------------------------------------------------
4383       real(kind=8) function RDIF()
4384
4385       use compare_data
4386       use geometry, only: dist
4387 !      implicit real*8 (a-h,o-z)
4388 !      include 'DIMENSIONS'
4389 !      include 'COMMON.CHAIN'
4390 !      include 'COMMON.DISTFIT'
4391       integer :: i,j,ind
4392       real(kind=8) :: suma,DIJ
4393 !      print *,'in rdif'
4394
4395       suma=0.0
4396       ind=0
4397       call chainbuild
4398       do i=1,nres-3
4399         do j=i+3,nres
4400           ind=ind+1
4401           if (w(ind).ne.0.0) then
4402             DIJ=DIST(i,j)
4403             suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind))
4404             DDD(ind)=DIJ
4405 !            print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma
4406           endif
4407         enddo
4408       enddo
4409
4410       RDIF=suma
4411       return
4412       end function RDIF
4413 !-----------------------------------------------------------------------------
4414       subroutine RDERIV
4415
4416       use compare_data
4417       use geometry_data
4418       use geometry, only:dist
4419 !      implicit real*8 (a-h,o-z)
4420 !      include 'DIMENSIONS'
4421 !      include 'COMMON.CHAIN'
4422 !      include 'COMMON.DISTFIT'
4423 !      include 'COMMON.GEO'
4424       integer :: i,j,k,l,I1,I2,IND
4425       real(kind=8),DIMENSION(3) :: E12,R13,R24,PRODU
4426
4427       DO I=1,NY
4428         DO J=1,NX
4429           DRDG(I,J)=0.0
4430         ENDDO
4431       ENDDO
4432       DO I=1,NX
4433         I1=I+1
4434         I2=I+2
4435         CALL VEC(I1,I2,E12)
4436         DO J=1,I
4437           DO K=1,3
4438             R13(K)=C(K,J)-C(K,I1)
4439           ENDDO
4440           DO K=I2+1,NRES
4441             DO L=1,3
4442               R24(L)=C(L,K)-C(L,I2)
4443             ENDDO
4444             IND=((J-1)*(2*NRES-J-6))/2+K-3
4445             PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2)
4446             PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3)
4447             PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1)
4448             DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K)
4449           ENDDO
4450         ENDDO
4451       ENDDO
4452       return
4453       end subroutine RDERIV
4454 !-----------------------------------------------------------------------------
4455       subroutine HEVAL
4456
4457       use compare_data
4458 !      implicit real*8 (a-h,o-z)
4459 !      include 'DIMENSIONS'
4460 !      include 'COMMON.CHAIN'
4461 !      include 'COMMON.DISTFIT'
4462       integer :: i,k,j
4463       real(kind=8) :: XI,HII,BKI,BKIWK,HIJ
4464
4465       DO I=1,NX
4466         XI=0.0
4467         HII=0.0
4468         DO K=1,NY
4469           BKI=DRDG(K,I)
4470           BKIWK=w(K)*BKI
4471           XI=XI+BKIWK*(D0(K)-DDD(K))
4472           HII=HII+BKI*BKIWK
4473         ENDDO
4474         H(I,I)=HII
4475         XX(I)=XI
4476         DO J=I+1,NX
4477           HIJ=0.0
4478           DO K=1,NY
4479             HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K)
4480           ENDDO
4481           H(I,J)=HIJ
4482           H(J,I)=HIJ
4483         ENDDO
4484       ENDDO
4485       return
4486       end subroutine HEVAL
4487 !-----------------------------------------------------------------------------
4488       subroutine VEC(I,J,U)
4489 !
4490       use geometry_data, only: C
4491 !  Find the unit vector from atom (I) to atom (J). Store in U.
4492 !
4493 !      implicit real*8 (a-h,o-z)
4494 !      include 'DIMENSIONS'
4495 !      include 'COMMON.CHAIN'
4496       integer :: I,J,K
4497       real(kind=8),DIMENSION(3) :: U
4498       real(kind=8) :: ANORM,UK
4499
4500       ANORM=0.0
4501       DO K=1,3
4502         UK=C(K,J)-C(K,I)
4503         ANORM=ANORM+UK*UK
4504         U(K)=UK
4505       ENDDO
4506       ANORM=SQRT(ANORM)
4507       DO K=1,3
4508         U(K)=U(K)/ANORM
4509       ENDDO
4510       return
4511       end subroutine VEC
4512 !-----------------------------------------------------------------------------
4513       subroutine TRANSFER(N,X1,X2)
4514
4515 !      implicit real*8 (a-h,o-z)
4516 !      include 'DIMENSIONS'
4517       integer :: N,I
4518       real(kind=8),DIMENSION(N) :: X1,X2
4519       DO 1 I=1,N
4520     1   X2(I)=X1(I)
4521       return
4522       end subroutine TRANSFER
4523 !-----------------------------------------------------------------------------
4524 !-----------------------------------------------------------------------------
4525       subroutine alloc_compare_arrays
4526
4527       maxres22=nres*(nres+1)/2
4528 ! common.dbase
4529 !      common /struct/ in io_common: read_threadbase
4530 !      allocate(cart_base !(3,maxres_base,maxseq)
4531 !      allocate(nres_base !(3,maxseq)
4532 !      allocate(str_nam !(maxseq)
4533 ! common.distfit
4534 !      COMMON /c_frag/ in io_conf: readpdb
4535       if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3)
4536       if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
4537 !      COMMON /WAGI/
4538       allocate(w(maxres22),d0(maxres22)) !(maxres22)
4539 !      COMMON /POCHODNE/
4540 !el      allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
4541       allocate(DDD(maxres22))   !(maxres22)
4542       allocate(H(nres,nres)) !(MAXRES,MAXRES)
4543       allocate(XX(nres)) !(MAXRES)
4544 !      COMMON /frozen/
4545       allocate(mask(nres)) !(maxres)
4546 ! common.thread
4547 !      common /thread/
4548       allocate(iexam(2,maxthread),ipatt(2,maxthread)) !(2,maxthread)
4549 !      common /thread1/
4550       allocate(ener0(n_ene+2,maxthread),ener(n_ene+2,maxthread)) !(n_ene+2,maxthread)
4551
4552       return
4553       end subroutine alloc_compare_arrays
4554 !-----------------------------------------------------------------------------
4555 #endif
4556 !-----------------------------------------------------------------------------
4557       end module compare