source/unres/src_MD-M/test.F

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