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