source/unres/compare.F90

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