source/unres/compare.F90

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