source/unres/compare.F90

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