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