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