source/unres/src_CSA/test.F

   1       subroutine test
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.GEO'
   8       include 'COMMON.VAR'
   9       include 'COMMON.INTERACT'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.DISTFIT'
  12       include 'COMMON.SBRIDGE'
  13       include 'COMMON.CONTROL'
  14       include 'COMMON.FFIELD'
  15       include 'COMMON.MINIM'
  16       include 'COMMON.CHAIN'
  17       double precision time0,time1
  18       double precision energy(0:n_ene),ee
  19       double precision var(maxvar),var1(maxvar)
  20       integer j1,j2
  21       logical debug,accepted
  22       debug=.true.
  23
  24
  25       call geom_to_var(nvar,var1)
  26       call chainbuild
  27       call etotal(energy(0))
  28       etot=energy(0)
  29       call rmsd(rms)
  30       write(iout,*) 'etot=',0,etot,rms
  31       call secondary2(.false.)
  32
  33       call write_pdb(0,'first structure',etot)
  34
  35       j1=13
  36       j2=21
  37       da=180.0*deg2rad
  38
  39
  40
  41        temp=3000.0d0
  42        betbol=1.0D0/(1.9858D-3*temp)
  43        jr=iran_num(j1,j2)
  44        d=ran_number(-pi,pi)
  45 c       phi(jr)=pinorm(phi(jr)+d)
  46        call chainbuild
  47        call etotal(energy(0))
  48        etot0=energy(0)
  49        call rmsd(rms)
  50        write(iout,*) 'etot=',1,etot0,rms
  51        call write_pdb(1,'perturb structure',etot0)
  52
  53       do i=2,500,2
  54        jr=iran_num(j1,j2)
  55        d=ran_number(-da,da)
  56        phiold=phi(jr)
  57        phi(jr)=pinorm(phi(jr)+d)
  58        call chainbuild
  59        call etotal(energy(0))
  60        etot=energy(0)
  61
  62        if (etot.lt.etot0) then
  63           accepted=.true.
  64        else
  65           accepted=.false.
  66           xxr=ran_number(0.0D0,1.0D0)
  67           xxh=betbol*(etot-etot0)
  68           if (xxh.lt.50.0D0) then
  69             xxh=dexp(-xxh)
  70             if (xxh.gt.xxr) accepted=.true.
  71           endif
  72        endif
  73        accepted=.true.
  74 c       print *,etot0,etot,accepted
  75        if (accepted) then
  76           etot0=etot
  77           call rmsd(rms)
  78           write(iout,*) 'etot=',i,etot,rms
  79           call write_pdb(i,'MC structure',etot)
  80 c minimize
  81 c        call geom_to_var(nvar,var1)
  82         call sc_move(2,nres-1,1,10d0,nft_sc,etot)
  83         call geom_to_var(nvar,var)
  84         call minimize(etot,var,iretcode,nfun)
  85         write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
  86         call var_to_geom(nvar,var)
  87         call chainbuild
  88         call rmsd(rms)
  89         write(iout,*) 'etot mcm=',i,etot,rms
  90         call write_pdb(i+1,'MCM structure',etot)
  91         call var_to_geom(nvar,var1)
  92 c --------
  93        else
  94           phi(jr)=phiold
  95        endif
  96       enddo
  97
  98 c minimize
  99 c       call sc_move(2,nres-1,1,10d0,nft_sc,etot)
 100 c       call geom_to_var(nvar,var)
 101 c
 102 c       call chainbuild
 103 c       call write_pdb(998 ,'sc min',etot)
 104 c
 105 c       call minimize(etot,var,iretcode,nfun)
 106 c       write(iout,*)'------------------------------------------------'
 107 c       write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
 108 c
 109 c       call var_to_geom(nvar,var)
 110 c       call chainbuild
 111 c       call write_pdb(999,'full min',etot)
 112
 113
 114       return
 115       end
 116
 117
 118       subroutine test_n16
 119       implicit real*8 (a-h,o-z)
 120       include 'DIMENSIONS'
 121 #ifdef MPI
 122       include 'mpif.h'
 123 #endif
 124       include 'COMMON.GEO'
 125       include 'COMMON.VAR'
 126       include 'COMMON.INTERACT'
 127       include 'COMMON.IOUNITS'
 128       include 'COMMON.DISTFIT'
 129       include 'COMMON.SBRIDGE'
 130       include 'COMMON.CONTROL'
 131       include 'COMMON.FFIELD'
 132       include 'COMMON.MINIM'
 133       include 'COMMON.CHAIN'
 134       double precision time0,time1
 135       double precision energy(0:n_ene),ee
 136       double precision var(maxvar),var1(maxvar)
 137       integer jdata(5)
 138       logical debug
 139       debug=.true.
 140
 141 c
 142       call geom_to_var(nvar,var1)
 143       call chainbuild
 144       call etotal(energy(0))
 145       etot=energy(0)
 146       write(iout,*) nnt,nct,etot
 147       call write_pdb(1,'first structure',etot)
 148       call secondary2(.true.)
 149
 150       do i=1,4
 151         jdata(i)=bfrag(i,2)
 152       enddo
 153
 154       DO ij=1,4
 155        ieval=0
 156        jdata(5)=ij
 157        call var_to_geom(nvar,var1)
 158        write(iout,*) 'N16 test',(jdata(i),i=1,5)
 159        call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5)
 160      &               ,ieval,ij)
 161        call geom_to_var(nvar,var)
 162
 163       if (minim) then
 164 #ifdef MPI
 165        time0=MPI_WTIME()
 166 #else
 167        time0=tcpu()
 168 #endif
 169        call minimize(etot,var,iretcode,nfun)
 170        write(iout,*)'------------------------------------------------'
 171        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
 172      &  '+ DIST eval',ieval
 173
 174 #ifdef MPI
 175        time1=MPI_WTIME()
 176 #else
 177        time1=tcpu()
 178 #endif
 179        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
 180      &         nfun/(time1-time0),' eval/s'
 181
 182        call var_to_geom(nvar,var)
 183        call chainbuild
 184        call write_pdb(ij*100+99,'full min',etot)
 185       endif
 186
 187
 188       ENDDO
 189
 190       return
 191       end
 192
 193
 194       subroutine test_local
 195       implicit real*8 (a-h,o-z)
 196       include 'DIMENSIONS'
 197       include 'COMMON.GEO'
 198       include 'COMMON.VAR'
 199       include 'COMMON.INTERACT'
 200       include 'COMMON.IOUNITS'
 201       double precision time0,time1
 202       double precision energy(0:n_ene),ee
 203       double precision varia(maxvar)
 204 c
 205       call chainbuild
 206 c      call geom_to_var(nvar,varia)
 207       call write_pdb(1,'first structure',0d0)
 208
 209       call etotal(energy(0))
 210       etot=energy(0)
 211       write(iout,*) nnt,nct,etot
 212
 213       write(iout,*) 'calling sc_move'
 214       call sc_move(nnt,nct,5,10d0,nft_sc,etot)
 215       write(iout,*) nft_sc,etot
 216       call write_pdb(2,'second structure',etot)
 217
 218       write(iout,*) 'calling local_move'
 219       call local_move_init(.false.)
 220       call local_move(24,29,20d0,50d0)
 221       call chainbuild
 222       call write_pdb(3,'third structure',etot)
 223
 224       write(iout,*) 'calling sc_move'
 225       call sc_move(24,29,5,10d0,nft_sc,etot)
 226       write(iout,*) nft_sc,etot
 227       call write_pdb(2,'last structure',etot)
 228
 229
 230       return
 231       end
 232
 233       subroutine test_sc
 234       implicit real*8 (a-h,o-z)
 235       include 'DIMENSIONS'
 236       include 'COMMON.GEO'
 237       include 'COMMON.VAR'
 238       include 'COMMON.INTERACT'
 239       include 'COMMON.IOUNITS'
 240       double precision time0,time1
 241       double precision energy(0:n_ene),ee
 242       double precision varia(maxvar)
 243 c
 244       call chainbuild
 245 c      call geom_to_var(nvar,varia)
 246       call write_pdb(1,'first structure',0d0)
 247
 248       call etotal(energy(0))
 249       etot=energy(0)
 250       write(iout,*) nnt,nct,etot
 251
 252       write(iout,*) 'calling sc_move'
 253
 254       call sc_move(nnt,nct,5,10d0,nft_sc,etot)
 255       write(iout,*) nft_sc,etot
 256       call write_pdb(2,'second structure',etot)
 257
 258       write(iout,*) 'calling sc_move 2nd time'
 259
 260       call sc_move(nnt,nct,5,1d0,nft_sc,etot)
 261       write(iout,*) nft_sc,etot
 262       call write_pdb(3,'last structure',etot)
 263       return
 264       end
 265 c--------------------------------------------------------
 266       subroutine bgrow(bstrand,nbstrand,in,ind,new)
 267       implicit real*8 (a-h,o-z)
 268       include 'DIMENSIONS'
 269       include 'COMMON.CHAIN'
 270       integer bstrand(maxres/3,6)
 271
 272       ishift=iabs(bstrand(in,ind+4)-new)
 273
 274       print *,'bgrow',bstrand(in,ind+4),new,ishift
 275
 276       bstrand(in,ind)=new
 277
 278       if(ind.eq.1)then
 279         bstrand(nbstrand,5)=bstrand(nbstrand,1)
 280         do i=1,nbstrand-1
 281           IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
 282           if (bstrand(i,5).lt.bstrand(i,6)) then
 283             bstrand(i,5)=bstrand(i,5)-ishift
 284           else
 285             bstrand(i,5)=bstrand(i,5)+ishift
 286           endif
 287           ENDIF
 288         enddo
 289       else
 290         bstrand(nbstrand,6)=bstrand(nbstrand,2)
 291         do i=1,nbstrand-1
 292           IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
 293           if (bstrand(i,6).lt.bstrand(i,5)) then
 294             bstrand(i,6)=bstrand(i,6)-ishift
 295           else
 296             bstrand(i,6)=bstrand(i,6)+ishift
 297           endif
 298           ENDIF
 299         enddo
 300       endif
 301
 302
 303       return
 304       end
 305
 306
 307 c------------------------------------------
 308       subroutine test11
 309       implicit real*8 (a-h,o-z)
 310       include 'DIMENSIONS'
 311 #ifdef MPI
 312       include 'mpif.h'
 313 #endif
 314       include 'COMMON.GEO'
 315       include 'COMMON.CHAIN'
 316       include 'COMMON.IOUNITS'
 317       include 'COMMON.VAR'
 318       include 'COMMON.CONTROL'
 319       include 'COMMON.SBRIDGE'
 320       include 'COMMON.FFIELD'
 321       include 'COMMON.MINIM'
 322 c
 323       include 'COMMON.DISTFIT'
 324       integer if(20,maxres),nif,ifa(20)
 325       integer ibc(0:maxres,0:maxres),istrand(20)
 326       integer ibd(maxres),ifb(10,2),nifb,lifb(10),lifb0
 327       integer itmp(20,maxres)
 328       double precision time0,time1
 329       double precision energy(0:n_ene),ee
 330       double precision varia(maxvar),vorg(maxvar)
 331 c
 332       logical debug,ltest,usedbfrag(maxres/3)
 333       character*50 linia
 334 c
 335       integer betasheet(maxres),ibetasheet(maxres),nbetasheet
 336       integer bstrand(maxres/3,6),nbstrand
 337
 338 c------------------------
 339
 340       debug=.true.
 341 c------------------------
 342       nbstrand=0
 343       nbetasheet=0
 344       do i=1,nres
 345         betasheet(i)=0
 346         ibetasheet(i)=0
 347       enddo
 348       call geom_to_var(nvar,vorg)
 349       call secondary2(debug)
 350
 351       if (nbfrag.le.1) return
 352
 353       do i=1,nbfrag
 354          usedbfrag(i)=.false.
 355       enddo
 356
 357
 358       nbetasheet=nbetasheet+1
 359       nbstrand=2
 360       bstrand(1,1)=bfrag(1,1)
 361       bstrand(1,2)=bfrag(2,1)
 362       bstrand(1,3)=nbetasheet
 363       bstrand(1,4)=1
 364       bstrand(1,5)=bfrag(1,1)
 365       bstrand(1,6)=bfrag(2,1)
 366       do i=bfrag(1,1),bfrag(2,1)
 367         betasheet(i)=nbetasheet
 368         ibetasheet(i)=1
 369       enddo
 370 c
 371       bstrand(2,1)=bfrag(3,1)
 372       bstrand(2,2)=bfrag(4,1)
 373       bstrand(2,3)=nbetasheet
 374       bstrand(2,5)=bfrag(3,1)
 375       bstrand(2,6)=bfrag(4,1)
 376
 377       if (bfrag(3,1).le.bfrag(4,1)) then
 378         bstrand(2,4)=2
 379         do i=bfrag(3,1),bfrag(4,1)
 380           betasheet(i)=nbetasheet
 381           ibetasheet(i)=2
 382         enddo
 383       else
 384         bstrand(2,4)=-2
 385         do i=bfrag(4,1),bfrag(3,1)
 386           betasheet(i)=nbetasheet
 387           ibetasheet(i)=2
 388         enddo
 389       endif
 390
 391       iused_nbfrag=1
 392
 393       do while (iused_nbfrag.ne.nbfrag)
 394
 395       do j=2,nbfrag
 396
 397         IF (.not.usedbfrag(j)) THEN
 398
 399         write (*,*) j,(bfrag(i,j),i=1,4)
 400         do jk=6,1,-1
 401          write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand)
 402         enddo
 403         write (*,*) '------------------'
 404
 405
 406         if (bfrag(3,j).le.bfrag(4,j)) then
 407          do i=bfrag(3,j),bfrag(4,j)
 408           if(betasheet(i).eq.nbetasheet) then
 409             in=ibetasheet(i)
 410             do k=bfrag(3,j),bfrag(4,j)
 411               betasheet(k)=nbetasheet
 412               ibetasheet(k)=in
 413             enddo
 414             nbstrand=nbstrand+1
 415             usedbfrag(j)=.true.
 416             iused_nbfrag=iused_nbfrag+1
 417             do k=bfrag(1,j),bfrag(2,j)
 418               betasheet(k)=nbetasheet
 419               ibetasheet(k)=nbstrand
 420             enddo
 421             if (bstrand(in,4).lt.0) then
 422               bstrand(nbstrand,1)=bfrag(2,j)
 423               bstrand(nbstrand,2)=bfrag(1,j)
 424               bstrand(nbstrand,3)=nbetasheet
 425               bstrand(nbstrand,4)=-nbstrand
 426               bstrand(nbstrand,5)=bstrand(nbstrand,1)
 427               bstrand(nbstrand,6)=bstrand(nbstrand,2)
 428               if(bstrand(in,1).lt.bfrag(4,j)) then
 429                  call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
 430               else
 431                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+
 432      &               (bstrand(in,5)-bfrag(4,j))
 433               endif
 434               if(bstrand(in,2).gt.bfrag(3,j)) then
 435                  call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
 436               else
 437                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-
 438      &               (-bstrand(in,6)+bfrag(3,j))
 439               endif
 440             else
 441               bstrand(nbstrand,1)=bfrag(1,j)
 442               bstrand(nbstrand,2)=bfrag(2,j)
 443               bstrand(nbstrand,3)=nbetasheet
 444               bstrand(nbstrand,4)=nbstrand
 445               bstrand(nbstrand,5)=bstrand(nbstrand,1)
 446               bstrand(nbstrand,6)=bstrand(nbstrand,2)
 447               if(bstrand(in,1).gt.bfrag(3,j)) then
 448                  call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
 449               else
 450                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-
 451      &                (-bstrand(in,5)+bfrag(3,j))
 452               endif
 453               if(bstrand(in,2).lt.bfrag(4,j)) then
 454                  call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
 455               else
 456                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+
 457      &              (bstrand(in,6)-bfrag(4,j))
 458               endif
 459             endif
 460             goto 11
 461           endif
 462           if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then
 463             in=ibetasheet(bfrag(1,j)+i-bfrag(3,j))
 464             do k=bfrag(1,j),bfrag(2,j)
 465               betasheet(k)=nbetasheet
 466               ibetasheet(k)=in
 467             enddo
 468             nbstrand=nbstrand+1
 469             usedbfrag(j)=.true.
 470             iused_nbfrag=iused_nbfrag+1
 471             do k=bfrag(3,1),bfrag(4,1)
 472               betasheet(k)=nbetasheet
 473               ibetasheet(k)=nbstrand
 474             enddo
 475             if (bstrand(in,4).lt.0) then
 476               bstrand(nbstrand,1)=bfrag(4,j)
 477               bstrand(nbstrand,2)=bfrag(3,j)
 478               bstrand(nbstrand,3)=nbetasheet
 479               bstrand(nbstrand,4)=-nbstrand
 480               bstrand(nbstrand,5)=bstrand(nbstrand,1)
 481               bstrand(nbstrand,6)=bstrand(nbstrand,2)
 482               if(bstrand(in,1).lt.bfrag(2,j)) then
 483                  call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
 484               else
 485                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+
 486      &              (bstrand(in,5)-bfrag(2,j))
 487               endif
 488               if(bstrand(in,2).gt.bfrag(1,j)) then
 489                  call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
 490               else
 491                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-
 492      &              (-bstrand(in,6)+bfrag(1,j))
 493               endif
 494             else
 495               bstrand(nbstrand,1)=bfrag(3,j)
 496               bstrand(nbstrand,2)=bfrag(4,j)
 497               bstrand(nbstrand,3)=nbetasheet
 498               bstrand(nbstrand,4)=nbstrand
 499               bstrand(nbstrand,5)=bstrand(nbstrand,1)
 500               bstrand(nbstrand,6)=bstrand(nbstrand,2)
 501               if(bstrand(in,1).gt.bfrag(1,j)) then
 502                  call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
 503               else
 504                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-
 505      &              (-bstrand(in,5)+bfrag(1,j))
 506               endif
 507               if(bstrand(in,2).lt.bfrag(2,j)) then
 508                  call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
 509               else
 510                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+
 511      &               (bstrand(in,6)-bfrag(2,j))
 512               endif
 513             endif
 514             goto 11
 515           endif
 516          enddo
 517         else
 518          do i=bfrag(4,j),bfrag(3,j)
 519           if(betasheet(i).eq.nbetasheet) then
 520             in=ibetasheet(i)
 521             do k=bfrag(4,j),bfrag(3,j)
 522               betasheet(k)=nbetasheet
 523               ibetasheet(k)=in
 524             enddo
 525             nbstrand=nbstrand+1
 526             usedbfrag(j)=.true.
 527             iused_nbfrag=iused_nbfrag+1
 528             do k=bfrag(1,j),bfrag(2,j)
 529               betasheet(k)=nbetasheet
 530               ibetasheet(k)=nbstrand
 531             enddo
 532             if (bstrand(in,4).lt.0) then
 533               bstrand(nbstrand,1)=bfrag(1,j)
 534               bstrand(nbstrand,2)=bfrag(2,j)
 535               bstrand(nbstrand,3)=nbetasheet
 536               bstrand(nbstrand,4)=nbstrand
 537               bstrand(nbstrand,5)=bstrand(nbstrand,1)
 538               bstrand(nbstrand,6)=bstrand(nbstrand,2)
 539               if(bstrand(in,1).lt.bfrag(3,j)) then
 540                  call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
 541               else
 542                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-
 543      &              (bstrand(in,5)-bfrag(3,j))
 544               endif
 545               if(bstrand(in,2).gt.bfrag(4,j)) then
 546                  call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
 547               else
 548                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+
 549      &              (-bstrand(in,6)+bfrag(4,j))
 550               endif
 551             else
 552               bstrand(nbstrand,1)=bfrag(2,j)
 553               bstrand(nbstrand,2)=bfrag(1,j)
 554               bstrand(nbstrand,3)=nbetasheet
 555               bstrand(nbstrand,4)=-nbstrand
 556               bstrand(nbstrand,5)=bstrand(nbstrand,1)
 557               bstrand(nbstrand,6)=bstrand(nbstrand,2)
 558               if(bstrand(in,1).gt.bfrag(4,j)) then
 559                  call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
 560               else
 561                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+
 562      &              (-bstrand(in,5)+bfrag(4,j))
 563               endif
 564               if(bstrand(in,2).lt.bfrag(3,j)) then
 565                  call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
 566               else
 567                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-
 568      &             (bstrand(in,6)-bfrag(3,j))
 569               endif
 570             endif
 571             goto 11
 572           endif
 573           if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then
 574             in=ibetasheet(bfrag(2,j)-i+bfrag(4,j))
 575             do k=bfrag(1,j),bfrag(2,j)
 576               betasheet(k)=nbetasheet
 577               ibetasheet(k)=in
 578             enddo
 579             nbstrand=nbstrand+1
 580             usedbfrag(j)=.true.
 581             iused_nbfrag=iused_nbfrag+1
 582             do k=bfrag(4,j),bfrag(3,j)
 583               betasheet(k)=nbetasheet
 584               ibetasheet(k)=nbstrand
 585             enddo
 586             if (bstrand(in,4).lt.0) then
 587               bstrand(nbstrand,1)=bfrag(4,j)
 588               bstrand(nbstrand,2)=bfrag(3,j)
 589               bstrand(nbstrand,3)=nbetasheet
 590               bstrand(nbstrand,4)=nbstrand
 591               bstrand(nbstrand,5)=bstrand(nbstrand,1)
 592               bstrand(nbstrand,6)=bstrand(nbstrand,2)
 593               if(bstrand(in,1).lt.bfrag(2,j)) then
 594                  call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
 595               else
 596                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-
 597      &             (bstrand(in,5)-bfrag(2,j))
 598               endif
 599               if(bstrand(in,2).gt.bfrag(1,j)) then
 600                  call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
 601               else
 602                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+
 603      &             (-bstrand(in,6)+bfrag(1,j))
 604               endif
 605             else
 606               bstrand(nbstrand,1)=bfrag(3,j)
 607               bstrand(nbstrand,2)=bfrag(4,j)
 608               bstrand(nbstrand,3)=nbetasheet
 609               bstrand(nbstrand,4)=-nbstrand
 610               bstrand(nbstrand,5)=bstrand(nbstrand,1)
 611               bstrand(nbstrand,6)=bstrand(nbstrand,2)
 612               if(bstrand(in,1).gt.bfrag(1,j)) then
 613                  call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
 614               else
 615                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+
 616      &               (-bstrand(in,5)+bfrag(1,j))
 617               endif
 618               if(bstrand(in,2).lt.bfrag(2,j)) then
 619                  call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
 620               else
 621                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-
 622      &                (bstrand(in,6)-bfrag(2,j))
 623               endif
 624             endif
 625             goto 11
 626           endif
 627          enddo
 628         endif
 629
 630
 631
 632         ENDIF
 633       enddo
 634
 635       j=2
 636       do while (usedbfrag(j))
 637         j=j+1
 638       enddo
 639
 640       nbstrand=nbstrand+1
 641       nbetasheet=nbetasheet+1
 642       bstrand(nbstrand,1)=bfrag(1,j)
 643       bstrand(nbstrand,2)=bfrag(2,j)
 644       bstrand(nbstrand,3)=nbetasheet
 645       bstrand(nbstrand,5)=bfrag(1,j)
 646       bstrand(nbstrand,6)=bfrag(2,j)
 647
 648       bstrand(nbstrand,4)=nbstrand
 649       do i=bfrag(1,j),bfrag(2,j)
 650           betasheet(i)=nbetasheet
 651           ibetasheet(i)=nbstrand
 652       enddo
 653 c
 654       nbstrand=nbstrand+1
 655       bstrand(nbstrand,1)=bfrag(3,j)
 656       bstrand(nbstrand,2)=bfrag(4,j)
 657       bstrand(nbstrand,3)=nbetasheet
 658       bstrand(nbstrand,5)=bfrag(3,j)
 659       bstrand(nbstrand,6)=bfrag(4,j)
 660
 661       if (bfrag(3,j).le.bfrag(4,j)) then
 662         bstrand(nbstrand,4)=nbstrand
 663         do i=bfrag(3,j),bfrag(4,j)
 664           betasheet(i)=nbetasheet
 665           ibetasheet(i)=nbstrand
 666         enddo
 667       else
 668         bstrand(nbstrand,4)=-nbstrand
 669         do i=bfrag(4,j),bfrag(3,j)
 670           betasheet(i)=nbetasheet
 671           ibetasheet(i)=nbstrand
 672         enddo
 673       endif
 674
 675       iused_nbfrag=iused_nbfrag+1
 676       usedbfrag(j)=.true.
 677
 678
 679   11  continue
 680         do jk=6,1,-1
 681          write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand)
 682         enddo
 683
 684
 685       enddo
 686
 687       do i=1,nres
 688        if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i)
 689       enddo
 690       write(*,*)
 691       do j=6,1,-1
 692         write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand)
 693       enddo
 694
 695 c------------------------
 696       nifb=0
 697       do i=1,nbstrand
 698         do j=i+1,nbstrand
 699            if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or.
 700      &          iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then
 701              nifb=nifb+1
 702              ifb(nifb,1)=bstrand(i,4)
 703              ifb(nifb,2)=bstrand(j,4)
 704            endif
 705         enddo
 706       enddo
 707
 708       write(*,*)
 709       do i=1,nifb
 710          write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2)
 711       enddo
 712
 713       do i=1,nbstrand
 714            ifa(i)=bstrand(i,4)
 715       enddo
 716       write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand)
 717
 718       nif=iabs(bstrand(1,6)-bstrand(1,5))+1
 719       do j=2,nbstrand
 720        if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif)
 721      &    nif=iabs(bstrand(j,6)-bstrand(j,5))+1
 722       enddo
 723
 724       write(*,*) nif
 725       do i=1,nif
 726         do j=1,nbstrand
 727           if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6))
 728           if (if(j,i).gt.0) then
 729             if(betasheet(if(j,i)).eq.0 .or.
 730      &          ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0
 731           else
 732             if(j,i)=0
 733           endif
 734         enddo
 735         write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand)
 736       enddo
 737
 738 c      read (inp,*) (ifa(i),i=1,4)
 739 c      do i=1,nres
 740 c       read (inp,*,err=20,end=20) (if(j,i),j=1,4)
 741 c      enddo
 742 c 20   nif=i-1
 743        stop
 744 c------------------------
 745
 746       isa=4
 747       is=2*isa-1
 748       iconf=0
 749 cccccccccccccccccccccccccccccccccc
 750       DO ig=1,is**isa-1
 751 cccccccccccccccccccccccccccccccccc
 752
 753          ii=ig
 754          do j=1,is
 755            istrand(is-j+1)=int(ii/is**(is-j))
 756            ii=ii-istrand(is-j+1)*is**(is-j)
 757          enddo
 758          ltest=.true.
 759          do k=1,isa
 760            istrand(k)=istrand(k)+1
 761            if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1
 762          enddo
 763          do k=1,isa
 764            do l=1,isa
 765             if(istrand(k).eq.istrand(l).and.k.ne.l.or.
 766      &          istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false.
 767            enddo
 768          enddo
 769
 770          lifb0=1
 771          do m=1,nifb
 772            lifb(m)=0
 773            do k=1,isa-1
 774             if(
 775      &         ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or.
 776      &         ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or.
 777      &       -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or.
 778      &       -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1))
 779      &         lifb(m)=1
 780            enddo
 781          lifb0=lifb0*lifb(m)
 782          enddo
 783
 784          if (mod(isa,2).eq.0) then
 785           do k=isa/2+1,isa
 786            if (istrand(k).eq.1) ltest=.false.
 787           enddo
 788          else
 789           do k=(isa+1)/2+1,isa
 790            if (istrand(k).eq.1) ltest=.false.
 791           enddo
 792          endif
 793
 794          IF (ltest.and.lifb0.eq.1) THEN
 795              iconf=iconf+1
 796
 797              call var_to_geom(nvar,vorg)
 798
 799              write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
 800              write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
 801              write (linia,'(10i3)') (istrand(k),k=1,isa)
 802
 803       do i=1,nres
 804         do j=1,nres
 805          ibc(i,j)=0
 806         enddo
 807       enddo
 808
 809
 810       do i=1,4
 811        if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then
 812         do j=1,nif
 813          itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j)
 814         enddo
 815        else
 816         do j=1,nif
 817         itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1)
 818         enddo
 819        endif
 820       enddo
 821
 822       do i=1,nif
 823         write(*,*) (itmp(j,i),j=1,4)
 824       enddo
 825
 826       do i=1,nif
 827 c       ifa(1),ifa(2),ifa(3),ifa(4)
 828 c       if(1,i),if(2,i),if(3,i),if(4,i)
 829         do k=1,isa-1
 830          ltest=.false.
 831          do m=1,nifb
 832            if(
 833      &         ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or.
 834      &         ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or.
 835      &       -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or.
 836      &       -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1))
 837      &   then
 838            ltest=.true.
 839            goto 110
 840          endif
 841          enddo
 842   110     continue
 843          if (ltest) then
 844           ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1
 845          else
 846           ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2
 847          endif
 848 c
 849         if (k.lt.3)
 850      &   ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3
 851         if (k.lt.2)
 852      &   ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4
 853         enddo
 854       enddo
 855 c------------------------
 856
 857 c
 858 c  freeze sec.elements
 859 c
 860        do i=1,nres
 861          mask(i)=1
 862          mask_phi(i)=1
 863          mask_theta(i)=1
 864          mask_side(i)=1
 865        enddo
 866
 867        do j=1,nbfrag
 868         do i=bfrag(1,j),bfrag(2,j)
 869          mask(i)=0
 870          mask_phi(i)=0
 871          mask_theta(i)=0
 872         enddo
 873         if (bfrag(3,j).le.bfrag(4,j)) then
 874          do i=bfrag(3,j),bfrag(4,j)
 875           mask(i)=0
 876           mask_phi(i)=0
 877           mask_theta(i)=0
 878          enddo
 879         else
 880          do i=bfrag(4,j),bfrag(3,j)
 881           mask(i)=0
 882           mask_phi(i)=0
 883           mask_theta(i)=0
 884          enddo
 885         endif
 886        enddo
 887        do j=1,nhfrag
 888         do i=hfrag(1,j),hfrag(2,j)
 889          mask(i)=0
 890          mask_phi(i)=0
 891          mask_theta(i)=0
 892         enddo
 893        enddo
 894        mask_r=.true.
 895
 896 c------------------------
 897 c      generate constrains
 898 c
 899        nhpb0=nhpb
 900        call chainbuild
 901        ind=0
 902        do i=1,nres-3
 903          do j=i+3,nres
 904           ind=ind+1
 905           if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then
 906             d0(ind)=DIST(i,j)
 907             w(ind)=10.0
 908             nhpb=nhpb+1
 909             ihpb(nhpb)=i
 910             jhpb(nhpb)=j
 911             forcon(nhpb)=10.0
 912             dhpb(nhpb)=d0(ind)
 913           else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then
 914             d0(ind)=5.0
 915             w(ind)=10.0
 916             nhpb=nhpb+1
 917             ihpb(nhpb)=i
 918             jhpb(nhpb)=j
 919             forcon(nhpb)=10.0
 920             dhpb(nhpb)=d0(ind)
 921           else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then
 922             d0(ind)=11.0
 923             w(ind)=10.0
 924             nhpb=nhpb+1
 925             ihpb(nhpb)=i
 926             jhpb(nhpb)=j
 927             forcon(nhpb)=10.0
 928             dhpb(nhpb)=d0(ind)
 929           else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then
 930             d0(ind)=16.0
 931             w(ind)=10.0
 932             nhpb=nhpb+1
 933             ihpb(nhpb)=i
 934             jhpb(nhpb)=j
 935             forcon(nhpb)=10.0
 936             dhpb(nhpb)=d0(ind)
 937           else if ( ibc(i,j).gt.0 ) then
 938             d0(ind)=DIST(i,ibc(i,j))
 939             w(ind)=10.0
 940             nhpb=nhpb+1
 941             ihpb(nhpb)=i
 942             jhpb(nhpb)=j
 943             forcon(nhpb)=10.0
 944             dhpb(nhpb)=d0(ind)
 945           else if ( ibc(j,i).gt.0 ) then
 946             d0(ind)=DIST(ibc(j,i),j)
 947             w(ind)=10.0
 948             nhpb=nhpb+1
 949             ihpb(nhpb)=i
 950             jhpb(nhpb)=j
 951             forcon(nhpb)=10.0
 952             dhpb(nhpb)=d0(ind)
 953           else
 954             w(ind)=0.0
 955           endif
 956           dd(ind)=d0(ind)
 957          enddo
 958        enddo
 959        call hpb_partition
 960 cd--------------------------
 961
 962       write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),
 963      &        ibc(jhpb(i),ihpb(i)),' --',
 964      &        ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
 965
 966 cd      nhpb=0
 967 cd      goto 901
 968 c
 969 c
 970       call contact_cp_min(varia,ifun,iconf,linia,debug)
 971       if (minim) then
 972 #ifdef MPI
 973        time0=MPI_WTIME()
 974 #else
 975        time0=tcpu()
 976 #endif
 977        call minimize(etot,varia,iretcode,nfun)
 978        write(iout,*)'------------------------------------------------'
 979        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
 980      &  '+ DIST eval',ifun
 981
 982 #ifdef MPI
 983        time1=MPI_WTIME()
 984 #else
 985        time1=tcpu()
 986 #endif
 987        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
 988      &         nfun/(time1-time0),' eval/s'
 989
 990        write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa)
 991        call var_to_geom(nvar,varia)
 992        call chainbuild
 993        call write_pdb(900+iconf,linia,etot)
 994       endif
 995
 996       call etotal(energy(0))
 997       etot=energy(0)
 998       call enerprint(energy(0))
 999 cd      call intout
1000 cd      call briefout(0,etot)
1001 cd      call secondary2(.true.)
1002
1003  901  CONTINUE
1004 ctest      return
1005 cccccccccccccccccccccccccccccccccccc
1006       ENDIF
1007       ENDDO
1008 cccccccccccccccccccccccccccccccccccc
1009
1010       return
1011   10  write (iout,'(a)') 'Error reading test structure.'
1012       return
1013       end
1014 c--------------------------------------------------------
1015
1016       subroutine test3
1017       implicit real*8 (a-h,o-z)
1018       include 'DIMENSIONS'
1019 #ifdef MPI
1020       include 'mpif.h'
1021 #endif
1022       include 'COMMON.GEO'
1023       include 'COMMON.CHAIN'
1024       include 'COMMON.IOUNITS'
1025       include 'COMMON.VAR'
1026       include 'COMMON.CONTROL'
1027       include 'COMMON.SBRIDGE'
1028       include 'COMMON.FFIELD'
1029       include 'COMMON.MINIM'
1030 c
1031       include 'COMMON.DISTFIT'
1032       integer if(3,maxres),nif
1033       integer ibc(maxres,maxres),istrand(20)
1034       integer ibd(maxres),ifb(10,2),nifb,lifb(10),lifb0
1035       double precision time0,time1
1036       double precision energy(0:n_ene),ee
1037       double precision varia(maxvar)
1038 c
1039       logical debug,ltest
1040       character*50 linia
1041 c
1042       do i=1,nres
1043        read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i)
1044       enddo
1045  20   nif=i-1
1046       write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),
1047      &                                         i=1,nif)
1048
1049
1050 c------------------------
1051       call secondary2(debug)
1052 c------------------------
1053       do i=1,nres
1054         do j=1,nres
1055          ibc(i,j)=0
1056         enddo
1057       enddo
1058
1059 c
1060 c  freeze sec.elements and store indexes for beta constrains
1061 c
1062        do i=1,nres
1063          mask(i)=1
1064          mask_phi(i)=1
1065          mask_theta(i)=1
1066          mask_side(i)=1
1067        enddo
1068
1069        do j=1,nbfrag
1070         do i=bfrag(1,j),bfrag(2,j)
1071          mask(i)=0
1072          mask_phi(i)=0
1073          mask_theta(i)=0
1074         enddo
1075         if (bfrag(3,j).le.bfrag(4,j)) then
1076          do i=bfrag(3,j),bfrag(4,j)
1077           mask(i)=0
1078           mask_phi(i)=0
1079           mask_theta(i)=0
1080           ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1
1081          enddo
1082         else
1083          do i=bfrag(4,j),bfrag(3,j)
1084           mask(i)=0
1085           mask_phi(i)=0
1086           mask_theta(i)=0
1087           ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1
1088          enddo
1089         endif
1090        enddo
1091        do j=1,nhfrag
1092         do i=hfrag(1,j),hfrag(2,j)
1093          mask(i)=0
1094          mask_phi(i)=0
1095          mask_theta(i)=0
1096         enddo
1097        enddo
1098        mask_r=.true.
1099
1100
1101 c ---------------- test --------------
1102        do i=1,nif
1103          if (ibc(if(1,i),if(2,i)).eq.-1) then
1104              ibc(if(1,i),if(2,i))=if(3,i)
1105              ibc(if(1,i),if(3,i))=if(2,i)
1106          else if (ibc(if(2,i),if(1,i)).eq.-1) then
1107              ibc(if(2,i),if(1,i))=0
1108              ibc(if(1,i),if(2,i))=if(3,i)
1109              ibc(if(1,i),if(3,i))=if(2,i)
1110          else
1111              ibc(if(1,i),if(2,i))=if(3,i)
1112              ibc(if(1,i),if(3,i))=if(2,i)
1113          endif
1114        enddo
1115
1116        do i=1,nres
1117         do j=1,nres
1118          if (ibc(i,j).ne.0)  write(*,'(3i5)') i,j,ibc(i,j)
1119         enddo
1120        enddo
1121 c------------------------
1122        call chainbuild
1123        ind=0
1124        do i=1,nres-3
1125          do j=i+3,nres
1126           ind=ind+1
1127           if ( ibc(i,j).eq.-1 ) then
1128             d0(ind)=DIST(i,j)
1129             w(ind)=10.0
1130             nhpb=nhpb+1
1131             ihpb(nhpb)=i
1132             jhpb(nhpb)=j
1133             forcon(nhpb)=10.0
1134             dhpb(nhpb)=d0(ind)
1135           else if ( ibc(i,j).gt.0 ) then
1136             d0(ind)=DIST(i,ibc(i,j))
1137             w(ind)=10.0
1138             nhpb=nhpb+1
1139             ihpb(nhpb)=i
1140             jhpb(nhpb)=j
1141             forcon(nhpb)=10.0
1142             dhpb(nhpb)=d0(ind)
1143           else if ( ibc(j,i).gt.0 ) then
1144             d0(ind)=DIST(ibc(j,i),j)
1145             w(ind)=10.0
1146             nhpb=nhpb+1
1147             ihpb(nhpb)=i
1148             jhpb(nhpb)=j
1149             forcon(nhpb)=10.0
1150             dhpb(nhpb)=d0(ind)
1151           else
1152             w(ind)=0.0
1153           endif
1154          enddo
1155        enddo
1156        call hpb_partition
1157
1158 cd--------------------------
1159       write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),
1160      &        ibc(jhpb(i),ihpb(i)),' --',
1161      &        ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
1162
1163
1164       linia='dist'
1165       debug=.true.
1166       in_pdb=7
1167 c
1168       call contact_cp_min(varia,ieval,in_pdb,linia,debug)
1169       if (minim) then
1170 #ifdef MPI
1171        time0=MPI_WTIME()
1172 #else
1173        time0=tcpu()
1174 #endif
1175        call minimize(etot,varia,iretcode,nfun)
1176        write(iout,*)'------------------------------------------------'
1177        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
1178      &  '+ DIST eval',ieval
1179
1180 #ifdef MPI
1181        time1=MPI_WTIME()
1182 #else
1183        time1=tcpu()
1184 #endif
1185        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
1186      &         nfun/(time1-time0),' eval/s'
1187
1188
1189        call var_to_geom(nvar,varia)
1190        call chainbuild
1191        call write_pdb(999,'full min',etot)
1192       endif
1193
1194       call etotal(energy(0))
1195       etot=energy(0)
1196       call enerprint(energy(0))
1197       call intout
1198       call briefout(0,etot)
1199       call secondary2(.true.)
1200
1201       return
1202   10  write (iout,'(a)') 'Error reading test structure.'
1203       return
1204       end
1205
1206
1207
1208
1209       subroutine test__
1210       implicit real*8 (a-h,o-z)
1211       include 'DIMENSIONS'
1212 #ifdef MPI
1213       include 'mpif.h'
1214 #endif
1215       include 'COMMON.GEO'
1216       include 'COMMON.CHAIN'
1217       include 'COMMON.IOUNITS'
1218       include 'COMMON.VAR'
1219       include 'COMMON.CONTROL'
1220       include 'COMMON.SBRIDGE'
1221       include 'COMMON.FFIELD'
1222       include 'COMMON.MINIM'
1223 c
1224       include 'COMMON.DISTFIT'
1225       integer if(2,2),ind
1226       integer iff(maxres)
1227       double precision time0,time1
1228       double precision energy(0:n_ene),ee
1229       double precision theta2(maxres),phi2(maxres),alph2(maxres),
1230      &              omeg2(maxres),
1231      &              theta1(maxres),phi1(maxres),alph1(maxres),
1232      &              omeg1(maxres)
1233       double precision varia(maxvar),varia2(maxvar)
1234 c
1235
1236
1237       read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2)
1238       write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2)
1239       read (inp,*,err=10,end=10) (theta2(i),i=3,nres)
1240       read (inp,*,err=10,end=10) (phi2(i),i=4,nres)
1241       read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1)
1242       read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1)
1243       do i=1,nres
1244         theta2(i)=deg2rad*theta2(i)
1245         phi2(i)=deg2rad*phi2(i)
1246         alph2(i)=deg2rad*alph2(i)
1247         omeg2(i)=deg2rad*omeg2(i)
1248       enddo
1249       do i=1,nres
1250         theta1(i)=theta(i)
1251         phi1(i)=phi(i)
1252         alph1(i)=alph(i)
1253         omeg1(i)=omeg(i)
1254       enddo
1255
1256       do i=1,nres
1257        mask(i)=1
1258       enddo
1259
1260
1261 c------------------------
1262       do i=1,nres
1263          iff(i)=0
1264       enddo
1265       do j=1,2
1266         do i=if(j,1),if(j,2)
1267           iff(i)=1
1268         enddo
1269       enddo
1270
1271       call chainbuild
1272       call geom_to_var(nvar,varia)
1273       call write_pdb(1,'first structure',0d0)
1274
1275         call secondary(.true.)
1276
1277         call secondary2(.true.)
1278
1279         do j=1,nbfrag
1280            if ( (bfrag(3,j).lt.bfrag(4,j) .or.
1281      &          bfrag(4,j)-bfrag(2,j).gt.4) .and.
1282      &          bfrag(2,j)-bfrag(1,j).gt.3 ) then
1283              nn=nn+1
1284
1285              if (bfrag(3,j).lt.bfrag(4,j)) then
1286                write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)')
1287      &                     "select",bfrag(1,j)-1,"-",bfrag(2,j)-1
1288      &                         ,",",bfrag(3,j)-1,"-",bfrag(4,j)-1
1289              else
1290                write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)')
1291      &                     "select",bfrag(1,j)-1,"-",bfrag(2,j)-1
1292      &                         ,",",bfrag(4,j)-1,"-",bfrag(3,j)-1
1293              endif
1294            endif
1295         enddo
1296
1297       do i=1,nres
1298         theta(i)=theta2(i)
1299         phi(i)=phi2(i)
1300         alph(i)=alph2(i)
1301         omeg(i)=omeg2(i)
1302       enddo
1303
1304       call chainbuild
1305       call geom_to_var(nvar,varia2)
1306       call write_pdb(2,'second structure',0d0)
1307
1308
1309
1310 c-------------------------------------------------------
1311
1312       ifun=-1
1313       call contact_cp(varia,varia2,iff,ifun,7)
1314       if (minim) then
1315 #ifdef MPI
1316        time0=MPI_WTIME()
1317 #else
1318        time0=tcpu()
1319 #endif
1320        call minimize(etot,varia,iretcode,nfun)
1321        write(iout,*)'------------------------------------------------'
1322        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
1323      &  '+ DIST eval',ifun
1324
1325 #ifdef MPI
1326        time1=MPI_WTIME()
1327 #else
1328        time1=tcpu()
1329 #endif
1330        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
1331      &         nfun/(time1-time0),' eval/s'
1332
1333
1334        call var_to_geom(nvar,varia)
1335        call chainbuild
1336        call write_pdb(999,'full min',etot)
1337       endif
1338
1339       call etotal(energy(0))
1340       etot=energy(0)
1341       call enerprint(energy(0))
1342       call intout
1343       call briefout(0,etot)
1344
1345       return
1346   10  write (iout,'(a)') 'Error reading test structure.'
1347       return
1348       end
1349
1350 c-------------------------------------------------
1351 c-------------------------------------------------
1352
1353       subroutine secondary(lprint)
1354       implicit real*8 (a-h,o-z)
1355       include 'DIMENSIONS'
1356       include 'COMMON.CHAIN'
1357       include 'COMMON.IOUNITS'
1358       include 'COMMON.DISTFIT'
1359
1360       integer ncont,icont(2,maxres*maxres/2),isec(maxres,3)
1361       logical lprint,not_done
1362       real dcont(maxres*maxres/2),d
1363       real rcomp /7.0/
1364       real rbeta /5.2/
1365       real ralfa /5.2/
1366       real r310 /6.6/
1367       double precision xpi(3),xpj(3)
1368
1369
1370
1371       call chainbuild
1372 cd      call write_pdb(99,'sec structure',0d0)
1373       ncont=0
1374       nbfrag=0
1375       nhfrag=0
1376       do i=1,nres
1377         isec(i,1)=0
1378         isec(i,2)=0
1379         isec(i,3)=0
1380       enddo
1381
1382       do i=2,nres-3
1383         do k=1,3
1384           xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
1385         enddo
1386         do j=i+2,nres
1387           do k=1,3
1388              xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
1389           enddo
1390 cd        d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
1391 cd     &         (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
1392 cd     &         (c(3,i)-c(3,j))*(c(3,i)-c(3,j))
1393 cd          print *,'CA',i,j,d
1394           d =  (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) +
1395      &         (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) +
1396      &         (xpi(3)-xpj(3))*(xpi(3)-xpj(3))
1397           if ( d.lt.rcomp*rcomp) then
1398             ncont=ncont+1
1399             icont(1,ncont)=i
1400             icont(2,ncont)=j
1401             dcont(ncont)=sqrt(d)
1402           endif
1403         enddo
1404       enddo
1405       if (lprint) then
1406         write (iout,*)
1407         write (iout,'(a)') '#PP contact map distances:'
1408         do i=1,ncont
1409           write (iout,'(3i4,f10.5)')
1410      &     i,icont(1,i),icont(2,i),dcont(i)
1411         enddo
1412       endif
1413
1414 c finding parallel beta
1415 cd      write (iout,*) '------- looking for parallel beta -----------'
1416       nbeta=0
1417       nstrand=0
1418       do i=1,ncont
1419         i1=icont(1,i)
1420         j1=icont(2,i)
1421         if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and.
1422      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
1423      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and.
1424      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and.
1425      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and.
1426      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
1427      &     ) then
1428           ii1=i1
1429           jj1=j1
1430 cd         write (iout,*) i1,j1,dcont(i)
1431           not_done=.true.
1432           do while (not_done)
1433            i1=i1+1
1434            j1=j1+1
1435             do j=1,ncont
1436               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)
1437      &              .and. dcont(j).le.rbeta .and.
1438      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
1439      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and.
1440      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and.
1441      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and.
1442      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
1443      &                            ) goto 5
1444             enddo
1445             not_done=.false.
1446   5         continue
1447 cd            write (iout,*) i1,j1,dcont(j),not_done
1448           enddo
1449           j1=j1-1
1450           i1=i1-1
1451           if (i1-ii1.gt.1) then
1452             ii1=max0(ii1-1,1)
1453             jj1=max0(jj1-1,1)
1454             nbeta=nbeta+1
1455             if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
1456
1457             nbfrag=nbfrag+1
1458             bfrag(1,nbfrag)=ii1
1459             bfrag(2,nbfrag)=i1
1460             bfrag(3,nbfrag)=jj1
1461             bfrag(4,nbfrag)=j1
1462
1463             do ij=ii1,i1
1464              isec(ij,1)=isec(ij,1)+1
1465              isec(ij,1+isec(ij,1))=nbeta
1466             enddo
1467             do ij=jj1,j1
1468              isec(ij,1)=isec(ij,1)+1
1469              isec(ij,1+isec(ij,1))=nbeta
1470             enddo
1471
1472            if(lprint) then
1473             nstrand=nstrand+1
1474             if (nbeta.le.9) then
1475               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
1476      &          "DefPropRes 'strand",nstrand,
1477      &          "' 'num = ",ii1-1,"..",i1-1,"'"
1478             else
1479               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
1480      &          "DefPropRes 'strand",nstrand,
1481      &          "' 'num = ",ii1-1,"..",i1-1,"'"
1482             endif
1483             nstrand=nstrand+1
1484             if (nbeta.le.9) then
1485               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
1486      &          "DefPropRes 'strand",nstrand,
1487      &          "' 'num = ",jj1-1,"..",j1-1,"'"
1488             else
1489               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
1490      &          "DefPropRes 'strand",nstrand,
1491      &          "' 'num = ",jj1-1,"..",j1-1,"'"
1492             endif
1493               write(12,'(a8,4i4)')
1494      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
1495            endif
1496           endif
1497         endif
1498       enddo
1499
1500 c finding antiparallel beta
1501 cd      write (iout,*) '--------- looking for antiparallel beta ---------'
1502
1503       do i=1,ncont
1504         i1=icont(1,i)
1505         j1=icont(2,i)
1506         if (dcont(i).le.rbeta.and.
1507      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
1508      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and.
1509      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and.
1510      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and.
1511      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
1512      &     ) then
1513           ii1=i1
1514           jj1=j1
1515 cd          write (iout,*) i1,j1,dcont(i)
1516
1517           not_done=.true.
1518           do while (not_done)
1519            i1=i1+1
1520            j1=j1-1
1521             do j=1,ncont
1522               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
1523      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
1524      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and.
1525      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and.
1526      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and.
1527      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
1528      &           .and. dcont(j).le.rbeta ) goto 6
1529             enddo
1530             not_done=.false.
1531   6         continue
1532 cd            write (iout,*) i1,j1,dcont(j),not_done
1533           enddo
1534           i1=i1-1
1535           j1=j1+1
1536           if (i1-ii1.gt.1) then
1537             if(lprint)write (iout,*)'antiparallel beta',
1538      &                   nbeta,ii1-1,i1,jj1,j1-1
1539
1540             nbfrag=nbfrag+1
1541             bfrag(1,nbfrag)=max0(ii1-1,1)
1542             bfrag(2,nbfrag)=i1
1543             bfrag(3,nbfrag)=jj1
1544             bfrag(4,nbfrag)=max0(j1-1,1)
1545
1546             nbeta=nbeta+1
1547             iii1=max0(ii1-1,1)
1548             do ij=iii1,i1
1549              isec(ij,1)=isec(ij,1)+1
1550              isec(ij,1+isec(ij,1))=nbeta
1551             enddo
1552             jjj1=max0(j1-1,1)
1553             do ij=jjj1,jj1
1554              isec(ij,1)=isec(ij,1)+1
1555              isec(ij,1+isec(ij,1))=nbeta
1556             enddo
1557
1558
1559            if (lprint) then
1560             nstrand=nstrand+1
1561             if (nstrand.le.9) then
1562               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
1563      &          "DefPropRes 'strand",nstrand,
1564      &          "' 'num = ",ii1-2,"..",i1-1,"'"
1565             else
1566               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
1567      &          "DefPropRes 'strand",nstrand,
1568      &          "' 'num = ",ii1-2,"..",i1-1,"'"
1569             endif
1570             nstrand=nstrand+1
1571             if (nstrand.le.9) then
1572               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
1573      &          "DefPropRes 'strand",nstrand,
1574      &          "' 'num = ",j1-2,"..",jj1-1,"'"
1575             else
1576               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
1577      &          "DefPropRes 'strand",nstrand,
1578      &          "' 'num = ",j1-2,"..",jj1-1,"'"
1579             endif
1580               write(12,'(a8,4i4)')
1581      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
1582            endif
1583           endif
1584         endif
1585       enddo
1586
1587       if (nstrand.gt.0.and.lprint) then
1588         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
1589         do i=2,nstrand
1590          if (i.le.9) then
1591           write(12,'(a9,i1,$)') " | strand",i
1592          else
1593           write(12,'(a9,i2,$)') " | strand",i
1594          endif
1595         enddo
1596         write(12,'(a1)') "'"
1597       endif
1598
1599
1600 c finding alpha or 310 helix
1601
1602       nhelix=0
1603       do i=1,ncont
1604         i1=icont(1,i)
1605         j1=icont(2,i)
1606         if (j1.eq.i1+3.and.dcont(i).le.r310
1607      &     .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
1608 cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
1609 cd          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
1610           ii1=i1
1611           jj1=j1
1612           if (isec(ii1,1).eq.0) then
1613             not_done=.true.
1614           else
1615             not_done=.false.
1616           endif
1617           do while (not_done)
1618             i1=i1+1
1619             j1=j1+1
1620             do j=1,ncont
1621               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
1622             enddo
1623             not_done=.false.
1624   10        continue
1625 cd            write (iout,*) i1,j1,not_done
1626           enddo
1627           j1=j1-1
1628           if (j1-ii1.gt.4) then
1629             nhelix=nhelix+1
1630 cd            write (iout,*)'helix',nhelix,ii1,j1
1631
1632             nhfrag=nhfrag+1
1633             hfrag(1,nhfrag)=ii1
1634             hfrag(2,nhfrag)=max0(j1-1,1)
1635
1636             do ij=ii1,j1
1637              isec(ij,1)=-1
1638             enddo
1639            if (lprint) then
1640             write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
1641             if (nhelix.le.9) then
1642               write(12,'(a17,i1,a9,i3,a2,i3,a1)')
1643      &          "DefPropRes 'helix",nhelix,
1644      &          "' 'num = ",ii1-1,"..",j1-2,"'"
1645             else
1646               write(12,'(a17,i2,a9,i3,a2,i3,a1)')
1647      &          "DefPropRes 'helix",nhelix,
1648      &          "' 'num = ",ii1-1,"..",j1-2,"'"
1649             endif
1650            endif
1651           endif
1652         endif
1653       enddo
1654
1655       if (nhelix.gt.0.and.lprint) then
1656         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
1657         do i=2,nhelix
1658          if (nhelix.le.9) then
1659           write(12,'(a8,i1,$)') " | helix",i
1660          else
1661           write(12,'(a8,i2,$)') " | helix",i
1662          endif
1663         enddo
1664         write(12,'(a1)') "'"
1665       endif
1666
1667       if (lprint) then
1668        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
1669        write(12,'(a20)') "XMacStand ribbon.mac"
1670       endif
1671
1672
1673       return
1674       end
1675 c----------------------------------------------------------------------------
1676
1677       subroutine write_pdb(npdb,titelloc,ee)
1678       implicit real*8 (a-h,o-z)
1679       include 'DIMENSIONS'
1680       include 'COMMON.IOUNITS'
1681       character*50 titelloc1
1682       character*(*) titelloc
1683       character*3 zahl
1684       character*5 liczba5
1685       double precision ee
1686       integer npdb,ilen
1687       external ilen
1688
1689       titelloc1=titelloc
1690       lenpre=ilen(prefix)
1691       if (npdb.lt.1000) then
1692        call numstr(npdb,zahl)
1693        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
1694       else
1695         if (npdb.lt.10000) then
1696          write(liczba5,'(i1,i4)') 0,npdb
1697         else
1698          write(liczba5,'(i5)') npdb
1699         endif
1700         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
1701       endif
1702       call pdbout(ee,titelloc1,ipdb)
1703       close(ipdb)
1704       return
1705       end
1706
1707 c-----------------------------------------------------------
1708       subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
1709       implicit real*8 (a-h,o-z)
1710       include 'DIMENSIONS'
1711 #ifdef MPI
1712       include 'mpif.h'
1713 #endif
1714       include 'COMMON.SBRIDGE'
1715       include 'COMMON.FFIELD'
1716       include 'COMMON.IOUNITS'
1717       include 'COMMON.DISTFIT'
1718       include 'COMMON.VAR'
1719       include 'COMMON.CHAIN'
1720       include 'COMMON.MINIM'
1721
1722       character*50 linia
1723       integer nf,ij(4)
1724       double precision var(maxvar),var2(maxvar)
1725       double precision time0,time1
1726       integer iff(maxres),ieval
1727       double precision theta1(maxres),phi1(maxres),alph1(maxres),
1728      &                 omeg1(maxres)
1729
1730
1731        call var_to_geom(nvar,var)
1732        call chainbuild
1733        nhpb0=nhpb
1734        ind=0
1735        do i=1,nres-3
1736          do j=i+3,nres
1737           ind=ind+1
1738           if ( iff(i).eq.1.and.iff(j).eq.1 ) then
1739             d0(ind)=DIST(i,j)
1740             w(ind)=10.0
1741             nhpb=nhpb+1
1742             ihpb(nhpb)=i
1743             jhpb(nhpb)=j
1744             forcon(nhpb)=10.0
1745             dhpb(nhpb)=d0(ind)
1746           else
1747             w(ind)=0.0
1748           endif
1749          enddo
1750        enddo
1751        call hpb_partition
1752
1753        do i=1,nres
1754         theta1(i)=theta(i)
1755         phi1(i)=phi(i)
1756         alph1(i)=alph(i)
1757         omeg1(i)=omeg(i)
1758        enddo
1759
1760        call var_to_geom(nvar,var2)
1761
1762        do i=1,nres
1763           if ( iff(i).eq.1 ) then
1764               theta(i)=theta1(i)
1765               phi(i)=phi1(i)
1766               alph(i)=alph1(i)
1767               omeg(i)=omeg1(i)
1768           endif
1769        enddo
1770
1771        call chainbuild
1772
1773        NX=NRES-3
1774        NY=((NRES-4)*(NRES-5))/2
1775        call distfit(.true.,200)
1776
1777
1778        ipot0=ipot
1779        maxmin0=maxmin
1780        maxfun0=maxfun
1781        wstrain0=wstrain
1782
1783        ipot=6
1784        maxmin=2000
1785        maxfun=5000
1786        call geom_to_var(nvar,var)
1787        call minimize(etot,var,iretcode,nfun)
1788        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
1789
1790         call var_to_geom(nvar,var)
1791         call chainbuild
1792
1793
1794         iwsk=0
1795         nf=0
1796         if (iff(1).eq.1) then
1797           iwsk=1
1798           nf=nf+1
1799           ij(nf)=0
1800         endif
1801         do i=2,nres
1802            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
1803              iwsk=1
1804              nf=nf+1
1805              ij(nf)=i
1806            endif
1807            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
1808              iwsk=0
1809              nf=nf+1
1810              ij(nf)=i-1
1811            endif
1812         enddo
1813         if (iff(nres).eq.1) then
1814           nf=nf+1
1815           ij(nf)=nres
1816         endif
1817
1818
1819 cd        write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
1820 cd     &                     "select",ij(1),"-",ij(2),
1821 cd     &                         ",",ij(3),"-",ij(4)
1822 cd        call write_pdb(in_pdb,linia,etot)
1823
1824
1825        ipot=ipot0
1826        maxmin=maxmin0
1827        maxfun=maxfun0
1828        call minimize(etot,var,iretcode,nfun)
1829        ieval=nfun
1830
1831        nhpb= nhpb0
1832        link_start=1
1833        link_end=nhpb
1834        wstrain=wstrain0
1835
1836        return
1837        end
1838 c-----------------------------------------------------------
1839       subroutine contact_cp(var,var2,iff,ieval,in_pdb)
1840       implicit real*8 (a-h,o-z)
1841       include 'DIMENSIONS'
1842       include 'COMMON.SBRIDGE'
1843       include 'COMMON.FFIELD'
1844       include 'COMMON.IOUNITS'
1845       include 'COMMON.DISTFIT'
1846       include 'COMMON.VAR'
1847       include 'COMMON.CHAIN'
1848       include 'COMMON.MINIM'
1849
1850       character*50 linia
1851       integer nf,ij(4)
1852       double precision energy(0:n_ene)
1853       double precision var(maxvar),var2(maxvar)
1854       double precision time0,time1
1855       integer iff(maxres),ieval
1856       double precision theta1(maxres),phi1(maxres),alph1(maxres),
1857      &                 omeg1(maxres)
1858       logical debug
1859
1860       debug=.false.
1861 c      debug=.true.
1862       if (ieval.eq.-1) debug=.true.
1863
1864
1865 c
1866 c store selected dist. constrains from 1st structure
1867 c
1868 #ifdef OSF
1869 c     Intercept NaNs in the coordinates
1870 c      write(iout,*) (var(i),i=1,nvar)
1871       x_sum=0.D0
1872       do i=1,nvar
1873         x_sum=x_sum+var(i)
1874       enddo
1875       if (x_sum.ne.x_sum) then
1876         write(iout,*)" *** contact_cp : Found NaN in coordinates"
1877         call flush(iout)
1878         print *," *** contact_cp : Found NaN in coordinates"
1879         return
1880       endif
1881 #endif
1882
1883
1884        call var_to_geom(nvar,var)
1885        call chainbuild
1886        nhpb0=nhpb
1887        ind=0
1888        do i=1,nres-3
1889          do j=i+3,nres
1890           ind=ind+1
1891           if ( iff(i).eq.1.and.iff(j).eq.1 ) then
1892             d0(ind)=DIST(i,j)
1893             w(ind)=10.0
1894             nhpb=nhpb+1
1895             ihpb(nhpb)=i
1896             jhpb(nhpb)=j
1897             forcon(nhpb)=10.0
1898             dhpb(nhpb)=d0(ind)
1899           else
1900             w(ind)=0.0
1901           endif
1902          enddo
1903        enddo
1904        call hpb_partition
1905
1906        do i=1,nres
1907         theta1(i)=theta(i)
1908         phi1(i)=phi(i)
1909         alph1(i)=alph(i)
1910         omeg1(i)=omeg(i)
1911        enddo
1912
1913 c
1914 c  freeze sec.elements from 2nd structure
1915 c
1916        do i=1,nres
1917          mask_phi(i)=1
1918          mask_theta(i)=1
1919          mask_side(i)=1
1920        enddo
1921
1922        call var_to_geom(nvar,var2)
1923        call secondary2(debug)
1924        do j=1,nbfrag
1925         do i=bfrag(1,j),bfrag(2,j)
1926          mask(i)=0
1927          mask_phi(i)=0
1928          mask_theta(i)=0
1929         enddo
1930         if (bfrag(3,j).le.bfrag(4,j)) then
1931          do i=bfrag(3,j),bfrag(4,j)
1932           mask(i)=0
1933           mask_phi(i)=0
1934           mask_theta(i)=0
1935          enddo
1936         else
1937          do i=bfrag(4,j),bfrag(3,j)
1938           mask(i)=0
1939           mask_phi(i)=0
1940           mask_theta(i)=0
1941          enddo
1942         endif
1943        enddo
1944        do j=1,nhfrag
1945         do i=hfrag(1,j),hfrag(2,j)
1946          mask(i)=0
1947          mask_phi(i)=0
1948          mask_theta(i)=0
1949         enddo
1950        enddo
1951        mask_r=.true.
1952
1953 c
1954 c      copy selected res from 1st to 2nd structure
1955 c
1956
1957        do i=1,nres
1958           if ( iff(i).eq.1 ) then
1959               theta(i)=theta1(i)
1960               phi(i)=phi1(i)
1961               alph(i)=alph1(i)
1962               omeg(i)=omeg1(i)
1963           endif
1964        enddo
1965
1966       if(debug) then
1967 c
1968 c     prepare description in linia variable
1969 c
1970         iwsk=0
1971         nf=0
1972         if (iff(1).eq.1) then
1973           iwsk=1
1974           nf=nf+1
1975           ij(nf)=1
1976         endif
1977         do i=2,nres
1978            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
1979              iwsk=1
1980              nf=nf+1
1981              ij(nf)=i
1982            endif
1983            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
1984              iwsk=0
1985              nf=nf+1
1986              ij(nf)=i-1
1987            endif
1988         enddo
1989         if (iff(nres).eq.1) then
1990           nf=nf+1
1991           ij(nf)=nres
1992         endif
1993
1994         write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
1995      &                     "SELECT",ij(1)-1,"-",ij(2)-1,
1996      &                         ",",ij(3)-1,"-",ij(4)-1
1997
1998       endif
1999 c
2000 c     run optimization
2001 c
2002       call contact_cp_min(var,ieval,in_pdb,linia,debug)
2003
2004       return
2005       end
2006
2007       subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
2008 c
2009 c    input : theta,phi,alph,omeg,in_pdb,linia,debug
2010 c    output : var,ieval
2011 c
2012       implicit real*8 (a-h,o-z)
2013       include 'DIMENSIONS'
2014 #ifdef MPI
2015       include 'mpif.h'
2016 #endif
2017       include 'COMMON.SBRIDGE'
2018       include 'COMMON.FFIELD'
2019       include 'COMMON.IOUNITS'
2020       include 'COMMON.DISTFIT'
2021       include 'COMMON.VAR'
2022       include 'COMMON.CHAIN'
2023       include 'COMMON.MINIM'
2024
2025       character*50 linia
2026       integer nf,ij(4)
2027       double precision energy(0:n_ene)
2028       double precision var(maxvar)
2029       double precision time0,time1
2030       integer ieval,info(3)
2031       logical debug,fail,check_var,reduce,change
2032
2033        write(iout,'(a20,i6,a20)')
2034      &             '------------------',in_pdb,'-------------------'
2035
2036        if (debug) then
2037         call chainbuild
2038         call write_pdb(1000+in_pdb,'combined structure',0d0)
2039 #ifdef MPI
2040        time0=MPI_WTIME()
2041 #else
2042        time0=tcpu()
2043 #endif
2044        endif
2045
2046 c
2047 c     run optimization of distances
2048 c
2049 c     uses d0(),w() and mask() for frozen 2D
2050 c
2051 ctest---------------------------------------------
2052 ctest       NX=NRES-3
2053 ctest       NY=((NRES-4)*(NRES-5))/2
2054 ctest       call distfit(debug,5000)
2055
2056        do i=1,nres
2057          mask_side(i)=0
2058        enddo
2059
2060        ipot01=ipot
2061        maxmin01=maxmin
2062        maxfun01=maxfun
2063 c       wstrain01=wstrain
2064        wsc01=wsc
2065        wscp01=wscp
2066        welec01=welec
2067        wvdwpp01=wvdwpp
2068 c      wang01=wang
2069        wscloc01=wscloc
2070        wtor01=wtor
2071        wtor_d01=wtor_d
2072
2073        ipot=6
2074        maxmin=2000
2075        maxfun=4000
2076 c       wstrain=1.0
2077        wsc=0.0
2078        wscp=0.0
2079        welec=0.0
2080        wvdwpp=0.0
2081 c      wang=0.0
2082        wscloc=0.0
2083        wtor=0.0
2084        wtor_d=0.0
2085
2086        call geom_to_var(nvar,var)
2087 cde       change=reduce(var)
2088        if (check_var(var,info)) then
2089           write(iout,*) 'cp_min error in input'
2090           print *,'cp_min error in input'
2091           return
2092        endif
2093
2094 cd       call etotal(energy(0))
2095 cd       call enerprint(energy(0))
2096 cd       call check_eint
2097
2098 #ifdef MPI
2099        time0=MPI_WTIME()
2100 #else
2101        time0=tcpu()
2102 #endif
2103 cdtest       call minimize(etot,var,iretcode,nfun)
2104 cdtest       write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
2105 #ifdef MPI
2106        time1=MPI_WTIME()
2107 #else
2108        time1=tcpu()
2109 #endif
2110
2111 cd       call etotal(energy(0))
2112 cd       call enerprint(energy(0))
2113 cd       call check_eint
2114
2115        do i=1,nres
2116          mask_side(i)=1
2117        enddo
2118
2119        ipot=ipot01
2120        maxmin=maxmin01
2121        maxfun=maxfun01
2122 c       wstrain=wstrain01
2123        wsc=wsc01
2124        wscp=wscp01
2125        welec=welec01
2126        wvdwpp=wvdwpp01
2127 c      wang=wang01
2128        wscloc=wscloc01
2129        wtor=wtor01
2130        wtor_d=wtor_d01
2131 ctest--------------------------------------------------
2132
2133        if(debug) then
2134 #ifdef MPI
2135        time1=MPI_WTIME()
2136 #else
2137        time1=tcpu()
2138 #endif
2139         write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
2140         call write_pdb(2000+in_pdb,'distfit structure',0d0)
2141        endif
2142
2143
2144        ipot0=ipot
2145        maxmin0=maxmin
2146        maxfun0=maxfun
2147        wstrain0=wstrain
2148 c
2149 c      run soft pot. optimization
2150 c         with constrains:
2151 c             nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
2152 c         and frozen 2D:
2153 c             mask_phi(),mask_theta(),mask_side(),mask_r
2154 c
2155        ipot=6
2156        maxmin=2000
2157        maxfun=4000
2158
2159 cde       change=reduce(var)
2160 cde       if (check_var(var,info)) write(iout,*) 'error before soft'
2161 #ifdef MPI
2162        time0=MPI_WTIME()
2163 #else
2164        time0=tcpu()
2165 #endif
2166        call minimize(etot,var,iretcode,nfun)
2167
2168        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2169 #ifdef MPI
2170        time1=MPI_WTIME()
2171 #else
2172        time1=tcpu()
2173 #endif
2174        write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2175      &         nfun/(time1-time0),' SOFT eval/s'
2176        if (debug) then
2177          call var_to_geom(nvar,var)
2178          call chainbuild
2179          call write_pdb(3000+in_pdb,'soft structure',etot)
2180        endif
2181 c
2182 c      run full UNRES optimization with constrains and frozen 2D
2183 c      the same variables as soft pot. optimizatio
2184 c
2185        ipot=ipot0
2186        maxmin=maxmin0
2187        maxfun=maxfun0
2188 c
2189 c check overlaps before calling full UNRES minim
2190 c
2191        call var_to_geom(nvar,var)
2192        call chainbuild
2193        call etotal(energy(0))
2194 #ifdef OSF
2195        write(iout,*) 'N7 ',energy(0)
2196        if (energy(0).ne.energy(0)) then
2197         write(iout,*) 'N7 error - gives NaN',energy(0)
2198        endif
2199 #endif
2200        ieval=1
2201        if (energy(1).eq.1.0d20) then
2202          write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
2203          call overlap_sc(fail)
2204          if(.not.fail) then
2205            call etotal(energy(0))
2206            ieval=ieval+1
2207            write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
2208          else
2209            mask_r=.false.
2210            nhpb= nhpb0
2211            link_start=1
2212            link_end=nhpb
2213            wstrain=wstrain0
2214            return
2215          endif
2216        endif
2217        call flush(iout)
2218 c
2219 cdte       time0=MPI_WTIME()
2220 cde       change=reduce(var)
2221 cde       if (check_var(var,info)) then
2222 cde         write(iout,*) 'error before mask dist'
2223 cde         call var_to_geom(nvar,var)
2224 cde         call chainbuild
2225 cde         call write_pdb(10000+in_pdb,'before mask dist',etot)
2226 cde       endif
2227 cdte       call minimize(etot,var,iretcode,nfun)
2228 cdte       write(iout,*)'SUMSL MASK DIST return code is',iretcode,
2229 cdte     &                          ' eval ',nfun
2230 cdte       ieval=ieval+nfun
2231 cdte
2232 cdte       time1=MPI_WTIME()
2233 cdte       write (iout,'(a,f6.2,f8.2,a)')
2234 cdte     &        ' Time for mask dist min.',time1-time0,
2235 cdte     &         nfun/(time1-time0),'  eval/s'
2236 cdte       call flush(iout)
2237        if (debug) then
2238          call var_to_geom(nvar,var)
2239          call chainbuild
2240          call write_pdb(4000+in_pdb,'mask dist',etot)
2241        endif
2242 c
2243 c      switch off freezing of 2D and
2244 c      run full UNRES optimization with constrains
2245 c
2246        mask_r=.false.
2247 #ifdef MPI
2248        time0=MPI_WTIME()
2249 #else
2250        time0=tcpu()
2251 #endif
2252 cde       change=reduce(var)
2253 cde       if (check_var(var,info)) then
2254 cde         write(iout,*) 'error before dist'
2255 cde         call var_to_geom(nvar,var)
2256 cde         call chainbuild
2257 cde         call write_pdb(11000+in_pdb,'before dist',etot)
2258 cde       endif
2259
2260        call minimize(etot,var,iretcode,nfun)
2261
2262 cde        change=reduce(var)
2263 cde        if (check_var(var,info)) then
2264 cde          write(iout,*) 'error after dist',ico
2265 cde          call var_to_geom(nvar,var)
2266 cde          call chainbuild
2267 cde          call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
2268 cde        endif
2269        write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2270        ieval=ieval+nfun
2271
2272 #ifdef MPI
2273        time1=MPI_WTIME()
2274 #else
2275        time1=tcpu()
2276 #endif
2277        write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,
2278      &         nfun/(time1-time0),'  eval/s'
2279 cde       call etotal(energy(0))
2280 cde       write(iout,*) 'N7 after dist',energy(0)
2281        call flush(iout)
2282
2283        if (debug) then
2284         call var_to_geom(nvar,var)
2285         call chainbuild
2286         call write_pdb(in_pdb,linia,etot)
2287        endif
2288 c
2289 c      reset constrains
2290 c
2291        nhpb= nhpb0
2292        link_start=1
2293        link_end=nhpb
2294        wstrain=wstrain0
2295
2296        return
2297        end
2298 c--------------------------------------------------------
2299       subroutine softreg
2300       implicit real*8 (a-h,o-z)
2301       include 'DIMENSIONS'
2302 #ifdef MPI
2303       include 'mpif.h'
2304 #endif
2305       include 'COMMON.GEO'
2306       include 'COMMON.CHAIN'
2307       include 'COMMON.IOUNITS'
2308       include 'COMMON.VAR'
2309       include 'COMMON.CONTROL'
2310       include 'COMMON.SBRIDGE'
2311       include 'COMMON.FFIELD'
2312       include 'COMMON.MINIM'
2313       include 'COMMON.INTERACT'
2314 c
2315       include 'COMMON.DISTFIT'
2316       integer iff(maxres)
2317       double precision time0,time1
2318       double precision energy(0:n_ene),ee
2319       double precision var(maxvar)
2320       integer ieval
2321 c
2322       logical debug,ltest,fail
2323       character*50 linia
2324 c
2325       linia='test'
2326       debug=.true.
2327       in_pdb=0
2328
2329
2330
2331 c------------------------
2332 c
2333 c  freeze sec.elements
2334 c
2335        do i=1,nres
2336          mask_phi(i)=1
2337          mask_theta(i)=1
2338          mask_side(i)=1
2339          iff(i)=0
2340        enddo
2341
2342        do j=1,nbfrag
2343         do i=bfrag(1,j),bfrag(2,j)
2344          mask_phi(i)=0
2345          mask_theta(i)=0
2346          iff(i)=1
2347         enddo
2348         if (bfrag(3,j).le.bfrag(4,j)) then
2349          do i=bfrag(3,j),bfrag(4,j)
2350           mask_phi(i)=0
2351           mask_theta(i)=0
2352           iff(i)=1
2353          enddo
2354         else
2355          do i=bfrag(4,j),bfrag(3,j)
2356           mask_phi(i)=0
2357           mask_theta(i)=0
2358           iff(i)=1
2359          enddo
2360         endif
2361        enddo
2362        do j=1,nhfrag
2363         do i=hfrag(1,j),hfrag(2,j)
2364          mask_phi(i)=0
2365          mask_theta(i)=0
2366          iff(i)=1
2367         enddo
2368        enddo
2369        mask_r=.true.
2370
2371
2372
2373        nhpb0=nhpb
2374 c
2375 c store dist. constrains
2376 c
2377        do i=1,nres-3
2378          do j=i+3,nres
2379            if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2380             nhpb=nhpb+1
2381             ihpb(nhpb)=i
2382             jhpb(nhpb)=j
2383             forcon(nhpb)=0.1
2384             dhpb(nhpb)=DIST(i,j)
2385            endif
2386          enddo
2387        enddo
2388        call hpb_partition
2389
2390        if (debug) then
2391         call chainbuild
2392         call write_pdb(100+in_pdb,'input reg. structure',0d0)
2393        endif
2394
2395
2396        ipot0=ipot
2397        maxmin0=maxmin
2398        maxfun0=maxfun
2399        wstrain0=wstrain
2400        wang0=wang
2401 c
2402 c      run soft pot. optimization
2403 c
2404        ipot=6
2405        wang=3.0
2406        maxmin=2000
2407        maxfun=4000
2408        call geom_to_var(nvar,var)
2409 #ifdef MPI
2410        time0=MPI_WTIME()
2411 #else
2412        time0=tcpu()
2413 #endif
2414        call minimize(etot,var,iretcode,nfun)
2415
2416        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2417 #ifdef MPI
2418        time1=MPI_WTIME()
2419 #else
2420        time1=tcpu()
2421 #endif
2422        write (iout,'(a,f6.2,f8.2,a)')'  Time for soft min.',time1-time0,
2423      &         nfun/(time1-time0),' SOFT eval/s'
2424        if (debug) then
2425          call var_to_geom(nvar,var)
2426          call chainbuild
2427          call write_pdb(300+in_pdb,'soft structure',etot)
2428        endif
2429 c
2430 c      run full UNRES optimization with constrains and frozen 2D
2431 c      the same variables as soft pot. optimizatio
2432 c
2433        ipot=ipot0
2434        wang=wang0
2435        maxmin=maxmin0
2436        maxfun=maxfun0
2437 #ifdef MPI
2438        time0=MPI_WTIME()
2439 #else
2440        time0=tcpu()
2441 #endif
2442        call minimize(etot,var,iretcode,nfun)
2443        write(iout,*)'SUMSL MASK DIST return code is',iretcode,
2444      &                          ' eval ',nfun
2445        ieval=nfun
2446
2447 #ifdef MPI
2448        time1=MPI_WTIME()
2449 #else
2450        time1=tcpu()
2451 #endif
2452        write (iout,'(a,f6.2,f8.2,a)')
2453      &        '  Time for mask dist min.',time1-time0,
2454      &         nfun/(time1-time0),'  eval/s'
2455        if (debug) then
2456          call var_to_geom(nvar,var)
2457          call chainbuild
2458          call write_pdb(400+in_pdb,'mask & dist',etot)
2459        endif
2460 c
2461 c      switch off constrains and
2462 c      run full UNRES optimization with frozen 2D
2463 c
2464
2465 c
2466 c      reset constrains
2467 c
2468        nhpb_c=nhpb
2469        nhpb=nhpb0
2470        link_start=1
2471        link_end=nhpb
2472        wstrain=wstrain0
2473
2474
2475 #ifdef MPI
2476        time0=MPI_WTIME()
2477 #else
2478        time0=tcpu()
2479 #endif
2480        call minimize(etot,var,iretcode,nfun)
2481        write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
2482        ieval=ieval+nfun
2483
2484 #ifdef MPI
2485        time1=MPI_WTIME()
2486 #else
2487        time1=tcpu()
2488 #endif
2489        write (iout,'(a,f6.2,f8.2,a)')'  Time for mask min.',time1-time0,
2490      &         nfun/(time1-time0),'  eval/s'
2491
2492
2493        if (debug) then
2494         call var_to_geom(nvar,var)
2495         call chainbuild
2496         call write_pdb(500+in_pdb,'mask 2d frozen',etot)
2497        endif
2498
2499        mask_r=.false.
2500
2501
2502 c
2503 c      run full UNRES optimization with constrains and NO frozen 2D
2504 c
2505
2506        nhpb=nhpb_c
2507        link_start=1
2508        link_end=nhpb
2509        maxfun=maxfun0/5
2510
2511        do ico=1,5
2512
2513        wstrain=wstrain0/ico
2514
2515 #ifdef MPI
2516        time0=MPI_WTIME()
2517 #else
2518        time0=tcpu()
2519 #endif
2520        call minimize(etot,var,iretcode,nfun)
2521        write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2522      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2523      &                          ' eval ',nfun
2524        ieval=nfun
2525
2526 #ifdef MPI
2527        time1=MPI_WTIME()
2528 #else
2529        time1=tcpu()
2530 #endif
2531        write (iout,'(a,f6.2,f8.2,a)')
2532      &        '  Time for dist min.',time1-time0,
2533      &         nfun/(time1-time0),'  eval/s'
2534        if (debug) then
2535          call var_to_geom(nvar,var)
2536          call chainbuild
2537          call write_pdb(600+in_pdb+ico,'dist cons',etot)
2538        endif
2539
2540        enddo
2541 c
2542        nhpb=nhpb0
2543        link_start=1
2544        link_end=nhpb
2545        wstrain=wstrain0
2546        maxfun=maxfun0
2547
2548
2549 c
2550       if (minim) then
2551 #ifdef MPI
2552        time0=MPI_WTIME()
2553 #else
2554        time0=tcpu()
2555 #endif
2556        call minimize(etot,var,iretcode,nfun)
2557        write(iout,*)'------------------------------------------------'
2558        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
2559      &  '+ DIST eval',ieval
2560
2561 #ifdef MPI
2562        time1=MPI_WTIME()
2563 #else
2564        time1=tcpu()
2565 #endif
2566        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
2567      &         nfun/(time1-time0),' eval/s'
2568
2569
2570        call var_to_geom(nvar,var)
2571        call chainbuild
2572        call write_pdb(999,'full min',etot)
2573       endif
2574
2575       return
2576       end
2577
2578
2579       subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
2580       implicit real*8 (a-h,o-z)
2581       include 'DIMENSIONS'
2582 #ifdef MPI
2583       include 'mpif.h'
2584 #endif
2585       include 'COMMON.GEO'
2586       include 'COMMON.VAR'
2587       include 'COMMON.INTERACT'
2588       include 'COMMON.IOUNITS'
2589       include 'COMMON.DISTFIT'
2590       include 'COMMON.SBRIDGE'
2591       include 'COMMON.CONTROL'
2592       include 'COMMON.FFIELD'
2593       include 'COMMON.MINIM'
2594       include 'COMMON.CHAIN'
2595       double precision time0,time1
2596       double precision energy(0:n_ene),ee
2597       double precision var(maxvar)
2598       integer jdata(5),isec(maxres)
2599 c
2600       jdata(1)=i1
2601       jdata(2)=i2
2602       jdata(3)=i3
2603       jdata(4)=i4
2604       jdata(5)=i5
2605
2606       call secondary2(.false.)
2607
2608       do i=1,nres
2609           isec(i)=0
2610       enddo
2611       do j=1,nbfrag
2612        do i=bfrag(1,j),bfrag(2,j)
2613           isec(i)=1
2614        enddo
2615        do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
2616           isec(i)=1
2617        enddo
2618       enddo
2619       do j=1,nhfrag
2620        do i=hfrag(1,j),hfrag(2,j)
2621           isec(i)=2
2622        enddo
2623       enddo
2624
2625 c
2626 c cut strands at the ends
2627 c
2628       if (jdata(2)-jdata(1).gt.3) then
2629        jdata(1)=jdata(1)+1
2630        jdata(2)=jdata(2)-1
2631        if (jdata(3).lt.jdata(4)) then
2632           jdata(3)=jdata(3)+1
2633           jdata(4)=jdata(4)-1
2634        else
2635           jdata(3)=jdata(3)-1
2636           jdata(4)=jdata(4)+1
2637        endif
2638       endif
2639
2640 cv      call chainbuild
2641 cv      call etotal(energy(0))
2642 cv      etot=energy(0)
2643 cv      write(iout,*) nnt,nct,etot
2644 cv      call write_pdb(ij*100,'first structure',etot)
2645 cv      write(iout,*) 'N16 test',(jdata(i),i=1,5)
2646
2647 c------------------------
2648 c      generate constrains
2649 c
2650        ishift=jdata(5)-2
2651        if(ishift.eq.0) ishift=-2
2652        nhpb0=nhpb
2653        call chainbuild
2654        do i=jdata(1),jdata(2)
2655         isec(i)=-1
2656         if(jdata(4).gt.jdata(3))then
2657          do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
2658             isec(j)=-1
2659 cd            print *,i,j,j+ishift
2660             nhpb=nhpb+1
2661             ihpb(nhpb)=i
2662             jhpb(nhpb)=j
2663             forcon(nhpb)=1000.0
2664             dhpb(nhpb)=DIST(i,j+ishift)
2665          enddo
2666         else
2667          do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
2668             isec(j)=-1
2669 cd            print *,i,j,j+ishift
2670             nhpb=nhpb+1
2671             ihpb(nhpb)=i
2672             jhpb(nhpb)=j
2673             forcon(nhpb)=1000.0
2674             dhpb(nhpb)=DIST(i,j+ishift)
2675          enddo
2676         endif
2677        enddo
2678
2679        do i=nnt,nct-2
2680          do j=i+2,nct
2681            if(isec(i).gt.0.or.isec(j).gt.0) then
2682 cd            print *,i,j
2683             nhpb=nhpb+1
2684             ihpb(nhpb)=i
2685             jhpb(nhpb)=j
2686             forcon(nhpb)=0.1
2687             dhpb(nhpb)=DIST(i,j)
2688            endif
2689          enddo
2690        enddo
2691
2692        call hpb_partition
2693
2694        call geom_to_var(nvar,var)
2695        maxfun0=maxfun
2696        wstrain0=wstrain
2697        maxfun=4000/5
2698
2699        do ico=1,5
2700
2701         wstrain=wstrain0/ico
2702
2703         call minimize(etot,var,iretcode,nfun)
2704         write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2705      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2706      &                          ' eval ',nfun
2707         ieval=ieval+nfun
2708
2709        enddo
2710 c
2711        nhpb=nhpb0
2712        call hpb_partition
2713        wstrain=wstrain0
2714        maxfun=maxfun0
2715 c
2716 cd      print *,etot
2717       wscloc0=wscloc
2718       wscloc=10.0
2719       call sc_move(nnt,nct,100,100d0,nft_sc,etot)
2720       wscloc=wscloc0
2721 cv      call chainbuild
2722 cv      call etotal(energy(0))
2723 cv      etot=energy(0)
2724 cv      call write_pdb(ij*100+10,'sc_move',etot)
2725 cd      call intout
2726 cd      print *,nft_sc,etot
2727
2728       return
2729       end
2730
2731       subroutine beta_zip(i1,i2,ieval,ij)
2732       implicit real*8 (a-h,o-z)
2733       include 'DIMENSIONS'
2734 #ifdef MPI
2735       include 'mpif.h'
2736 #endif
2737       include 'COMMON.GEO'
2738       include 'COMMON.VAR'
2739       include 'COMMON.INTERACT'
2740       include 'COMMON.IOUNITS'
2741       include 'COMMON.DISTFIT'
2742       include 'COMMON.SBRIDGE'
2743       include 'COMMON.CONTROL'
2744       include 'COMMON.FFIELD'
2745       include 'COMMON.MINIM'
2746       include 'COMMON.CHAIN'
2747       double precision time0,time1
2748       double precision energy(0:n_ene),ee
2749       double precision var(maxvar)
2750       character*10 test
2751
2752 cv      call chainbuild
2753 cv      call etotal(energy(0))
2754 cv      etot=energy(0)
2755 cv      write(test,'(2i5)') i1,i2
2756 cv      call write_pdb(ij*100,test,etot)
2757 cv      write(iout,*) 'N17 test',i1,i2,etot,ij
2758
2759 c
2760 c      generate constrains
2761 c
2762        nhpb0=nhpb
2763        nhpb=nhpb+1
2764        ihpb(nhpb)=i1
2765        jhpb(nhpb)=i2
2766        forcon(nhpb)=1000.0
2767        dhpb(nhpb)=4.0
2768
2769        call hpb_partition
2770
2771        call geom_to_var(nvar,var)
2772        maxfun0=maxfun
2773        wstrain0=wstrain
2774        maxfun=1000/5
2775
2776        do ico=1,5
2777         wstrain=wstrain0/ico
2778         call minimize(etot,var,iretcode,nfun)
2779         write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2780      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2781      &                          ' eval ',nfun
2782         ieval=ieval+nfun
2783 c do not comment the next line
2784          call var_to_geom(nvar,var)
2785 cv         call chainbuild
2786 cv         call write_pdb(ij*100+ico,'dist cons',etot)
2787        enddo
2788
2789        nhpb=nhpb0
2790        call hpb_partition
2791        wstrain=wstrain0
2792        maxfun=maxfun0
2793
2794 cv      call etotal(energy(0))
2795 cv      etot=energy(0)
2796 cv      write(iout,*) 'N17 test end',i1,i2,etot,ij
2797
2798
2799       return
2800       end