source/unres/src_MD-M/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 cd       call write_pdb(3,'combined structure',0d0)
1773 cd       time0=MPI_WTIME()
1774
1775        NX=NRES-3
1776        NY=((NRES-4)*(NRES-5))/2
1777        call distfit(.true.,200)
1778
1779 cd       time1=MPI_WTIME()
1780 cd       write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
1781
1782        ipot0=ipot
1783        maxmin0=maxmin
1784        maxfun0=maxfun
1785        wstrain0=wstrain
1786
1787        ipot=6
1788        maxmin=2000
1789        maxfun=5000
1790        call geom_to_var(nvar,var)
1791 cd       time0=MPI_WTIME()
1792        call minimize(etot,var,iretcode,nfun)
1793        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
1794
1795 cd       time1=MPI_WTIME()
1796 cd       write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
1797 cd     &         nfun/(time1-time0),' SOFT eval/s'
1798         call var_to_geom(nvar,var)
1799         call chainbuild
1800
1801
1802         iwsk=0
1803         nf=0
1804         if (iff(1).eq.1) then
1805           iwsk=1
1806           nf=nf+1
1807           ij(nf)=0
1808         endif
1809         do i=2,nres
1810            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
1811              iwsk=1
1812              nf=nf+1
1813              ij(nf)=i
1814            endif
1815            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
1816              iwsk=0
1817              nf=nf+1
1818              ij(nf)=i-1
1819            endif
1820         enddo
1821         if (iff(nres).eq.1) then
1822           nf=nf+1
1823           ij(nf)=nres
1824         endif
1825
1826
1827 cd        write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
1828 cd     &                     "select",ij(1),"-",ij(2),
1829 cd     &                         ",",ij(3),"-",ij(4)
1830 cd        call write_pdb(in_pdb,linia,etot)
1831
1832
1833        ipot=ipot0
1834        maxmin=maxmin0
1835        maxfun=maxfun0
1836 cd       time0=MPI_WTIME()
1837        call minimize(etot,var,iretcode,nfun)
1838 cd       write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
1839        ieval=nfun
1840
1841 cd       time1=MPI_WTIME()
1842 cd       write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
1843 cd     &         nfun/(time1-time0),' eval/s'
1844 cd       call var_to_geom(nvar,var)
1845 cd       call chainbuild
1846 cd       call write_pdb(6,'dist structure',etot)
1847
1848
1849        nhpb= nhpb0
1850        link_start=1
1851        link_end=nhpb
1852        wstrain=wstrain0
1853
1854        return
1855        end
1856 c-----------------------------------------------------------
1857       subroutine contact_cp(var,var2,iff,ieval,in_pdb)
1858       implicit real*8 (a-h,o-z)
1859       include 'DIMENSIONS'
1860       include 'COMMON.SBRIDGE'
1861       include 'COMMON.FFIELD'
1862       include 'COMMON.IOUNITS'
1863       include 'COMMON.DISTFIT'
1864       include 'COMMON.VAR'
1865       include 'COMMON.CHAIN'
1866       include 'COMMON.MINIM'
1867
1868       character*50 linia
1869       integer nf,ij(4)
1870       double precision energy(0:n_ene)
1871       double precision var(maxvar),var2(maxvar)
1872       double precision time0,time1
1873       integer iff(maxres),ieval
1874       double precision theta1(maxres),phi1(maxres),alph1(maxres),
1875      &                 omeg1(maxres)
1876       logical debug
1877
1878       debug=.false.
1879 c      debug=.true.
1880       if (ieval.eq.-1) debug=.true.
1881
1882
1883 c
1884 c store selected dist. constrains from 1st structure
1885 c
1886 #ifdef OSF
1887 c     Intercept NaNs in the coordinates
1888 c      write(iout,*) (var(i),i=1,nvar)
1889       x_sum=0.D0
1890       do i=1,nvar
1891         x_sum=x_sum+var(i)
1892       enddo
1893       if (x_sum.ne.x_sum) then
1894         write(iout,*)" *** contact_cp : Found NaN in coordinates"
1895         call flush(iout)
1896         print *," *** contact_cp : Found NaN in coordinates"
1897         return
1898       endif
1899 #endif
1900
1901
1902        call var_to_geom(nvar,var)
1903        call chainbuild
1904        nhpb0=nhpb
1905        ind=0
1906        do i=1,nres-3
1907          do j=i+3,nres
1908           ind=ind+1
1909           if ( iff(i).eq.1.and.iff(j).eq.1 ) then
1910             d0(ind)=DIST(i,j)
1911             w(ind)=10.0
1912             nhpb=nhpb+1
1913             ihpb(nhpb)=i
1914             jhpb(nhpb)=j
1915             forcon(nhpb)=10.0
1916             dhpb(nhpb)=d0(ind)
1917           else
1918             w(ind)=0.0
1919           endif
1920          enddo
1921        enddo
1922        call hpb_partition
1923
1924        do i=1,nres
1925         theta1(i)=theta(i)
1926         phi1(i)=phi(i)
1927         alph1(i)=alph(i)
1928         omeg1(i)=omeg(i)
1929        enddo
1930
1931 c
1932 c  freeze sec.elements from 2nd structure
1933 c
1934        do i=1,nres
1935          mask_phi(i)=1
1936          mask_theta(i)=1
1937          mask_side(i)=1
1938        enddo
1939
1940        call var_to_geom(nvar,var2)
1941        call secondary2(debug)
1942        do j=1,nbfrag
1943         do i=bfrag(1,j),bfrag(2,j)
1944          mask(i)=0
1945          mask_phi(i)=0
1946          mask_theta(i)=0
1947         enddo
1948         if (bfrag(3,j).le.bfrag(4,j)) then
1949          do i=bfrag(3,j),bfrag(4,j)
1950           mask(i)=0
1951           mask_phi(i)=0
1952           mask_theta(i)=0
1953          enddo
1954         else
1955          do i=bfrag(4,j),bfrag(3,j)
1956           mask(i)=0
1957           mask_phi(i)=0
1958           mask_theta(i)=0
1959          enddo
1960         endif
1961        enddo
1962        do j=1,nhfrag
1963         do i=hfrag(1,j),hfrag(2,j)
1964          mask(i)=0
1965          mask_phi(i)=0
1966          mask_theta(i)=0
1967         enddo
1968        enddo
1969        mask_r=.true.
1970
1971 c
1972 c      copy selected res from 1st to 2nd structure
1973 c
1974
1975        do i=1,nres
1976           if ( iff(i).eq.1 ) then
1977               theta(i)=theta1(i)
1978               phi(i)=phi1(i)
1979               alph(i)=alph1(i)
1980               omeg(i)=omeg1(i)
1981           endif
1982        enddo
1983
1984       if(debug) then
1985 c
1986 c     prepare description in linia variable
1987 c
1988         iwsk=0
1989         nf=0
1990         if (iff(1).eq.1) then
1991           iwsk=1
1992           nf=nf+1
1993           ij(nf)=1
1994         endif
1995         do i=2,nres
1996            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
1997              iwsk=1
1998              nf=nf+1
1999              ij(nf)=i
2000            endif
2001            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2002              iwsk=0
2003              nf=nf+1
2004              ij(nf)=i-1
2005            endif
2006         enddo
2007         if (iff(nres).eq.1) then
2008           nf=nf+1
2009           ij(nf)=nres
2010         endif
2011
2012         write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
2013      &                     "SELECT",ij(1)-1,"-",ij(2)-1,
2014      &                         ",",ij(3)-1,"-",ij(4)-1
2015
2016       endif
2017 c
2018 c     run optimization
2019 c
2020       call contact_cp_min(var,ieval,in_pdb,linia,debug)
2021
2022       return
2023       end
2024
2025       subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
2026 c
2027 c    input : theta,phi,alph,omeg,in_pdb,linia,debug
2028 c    output : var,ieval
2029 c
2030       implicit real*8 (a-h,o-z)
2031       include 'DIMENSIONS'
2032 #ifdef MPI
2033       include 'mpif.h'
2034 #endif
2035       include 'COMMON.SBRIDGE'
2036       include 'COMMON.FFIELD'
2037       include 'COMMON.IOUNITS'
2038       include 'COMMON.DISTFIT'
2039       include 'COMMON.VAR'
2040       include 'COMMON.CHAIN'
2041       include 'COMMON.MINIM'
2042
2043       character*50 linia
2044       integer nf,ij(4)
2045       double precision energy(0:n_ene)
2046       double precision var(maxvar)
2047       double precision time0,time1
2048       integer ieval,info(3)
2049       logical debug,fail,check_var,reduce,change
2050
2051        write(iout,'(a20,i6,a20)')
2052      &             '------------------',in_pdb,'-------------------'
2053
2054        if (debug) then
2055         call chainbuild
2056         call write_pdb(1000+in_pdb,'combined structure',0d0)
2057 #ifdef MPI
2058         time0=MPI_WTIME()
2059 #else
2060         time0=tcpu()
2061 #endif
2062        endif
2063
2064 c
2065 c     run optimization of distances
2066 c
2067 c     uses d0(),w() and mask() for frozen 2D
2068 c
2069 ctest---------------------------------------------
2070 ctest       NX=NRES-3
2071 ctest       NY=((NRES-4)*(NRES-5))/2
2072 ctest       call distfit(debug,5000)
2073
2074        do i=1,nres
2075          mask_side(i)=0
2076        enddo
2077
2078        ipot01=ipot
2079        maxmin01=maxmin
2080        maxfun01=maxfun
2081 c       wstrain01=wstrain
2082        wsc01=wsc
2083        wscp01=wscp
2084        welec01=welec
2085        wvdwpp01=wvdwpp
2086 c      wang01=wang
2087        wscloc01=wscloc
2088        wtor01=wtor
2089        wtor_d01=wtor_d
2090
2091        ipot=6
2092        maxmin=2000
2093        maxfun=4000
2094 c       wstrain=1.0
2095        wsc=0.0
2096        wscp=0.0
2097        welec=0.0
2098        wvdwpp=0.0
2099 c      wang=0.0
2100        wscloc=0.0
2101        wtor=0.0
2102        wtor_d=0.0
2103
2104        call geom_to_var(nvar,var)
2105 cde       change=reduce(var)
2106        if (check_var(var,info)) then
2107           write(iout,*) 'cp_min error in input'
2108           print *,'cp_min error in input'
2109           return
2110        endif
2111
2112 cd       call etotal(energy(0))
2113 cd       call enerprint(energy(0))
2114 cd       call check_eint
2115
2116 #ifdef MPI
2117        time0=MPI_WTIME()
2118 #else
2119        time0=tcpu()
2120 #endif
2121 cdtest       call minimize(etot,var,iretcode,nfun)
2122 cdtest       write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
2123 #ifdef MPI
2124        time1=MPI_WTIME()
2125 #else
2126        time1=tcpu()
2127 #endif
2128
2129 cd       call etotal(energy(0))
2130 cd       call enerprint(energy(0))
2131 cd       call check_eint
2132
2133        do i=1,nres
2134          mask_side(i)=1
2135        enddo
2136
2137        ipot=ipot01
2138        maxmin=maxmin01
2139        maxfun=maxfun01
2140 c       wstrain=wstrain01
2141        wsc=wsc01
2142        wscp=wscp01
2143        welec=welec01
2144        wvdwpp=wvdwpp01
2145 c      wang=wang01
2146        wscloc=wscloc01
2147        wtor=wtor01
2148        wtor_d=wtor_d01
2149 ctest--------------------------------------------------
2150
2151        if(debug) then
2152 #ifdef MPI
2153         time1=MPI_WTIME()
2154 #else
2155         time1=tcpu()
2156 #endif
2157         write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
2158         call write_pdb(2000+in_pdb,'distfit structure',0d0)
2159        endif
2160
2161
2162        ipot0=ipot
2163        maxmin0=maxmin
2164        maxfun0=maxfun
2165        wstrain0=wstrain
2166 c
2167 c      run soft pot. optimization
2168 c         with constrains:
2169 c             nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
2170 c         and frozen 2D:
2171 c             mask_phi(),mask_theta(),mask_side(),mask_r
2172 c
2173        ipot=6
2174        maxmin=2000
2175        maxfun=4000
2176
2177 cde       change=reduce(var)
2178 cde       if (check_var(var,info)) write(iout,*) 'error before soft'
2179 #ifdef MPI
2180        time0=MPI_WTIME()
2181 #else
2182        time0=tcpu()
2183 #endif
2184        call minimize(etot,var,iretcode,nfun)
2185
2186        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2187 #ifdef MPI
2188        time1=MPI_WTIME()
2189 #else
2190        time1=tcpu()
2191 #endif
2192        write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2193      &         nfun/(time1-time0),' SOFT eval/s'
2194        if (debug) then
2195          call var_to_geom(nvar,var)
2196          call chainbuild
2197          call write_pdb(3000+in_pdb,'soft structure',etot)
2198        endif
2199 c
2200 c      run full UNRES optimization with constrains and frozen 2D
2201 c      the same variables as soft pot. optimizatio
2202 c
2203        ipot=ipot0
2204        maxmin=maxmin0
2205        maxfun=maxfun0
2206 c
2207 c check overlaps before calling full UNRES minim
2208 c
2209        call var_to_geom(nvar,var)
2210        call chainbuild
2211        call etotal(energy(0))
2212 #ifdef OSF
2213        write(iout,*) 'N7 ',energy(0)
2214        if (energy(0).ne.energy(0)) then
2215         write(iout,*) 'N7 error - gives NaN',energy(0)
2216        endif
2217 #endif
2218        ieval=1
2219        if (energy(1).eq.1.0d20) then
2220          write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
2221          call overlap_sc(fail)
2222          if(.not.fail) then
2223            call etotal(energy(0))
2224            ieval=ieval+1
2225            write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
2226          else
2227            mask_r=.false.
2228            nhpb= nhpb0
2229            link_start=1
2230            link_end=nhpb
2231            wstrain=wstrain0
2232            return
2233          endif
2234        endif
2235        call flush(iout)
2236 c
2237 cdte       time0=MPI_WTIME()
2238 cde       change=reduce(var)
2239 cde       if (check_var(var,info)) then
2240 cde         write(iout,*) 'error before mask dist'
2241 cde         call var_to_geom(nvar,var)
2242 cde         call chainbuild
2243 cde         call write_pdb(10000+in_pdb,'before mask dist',etot)
2244 cde       endif
2245 cdte       call minimize(etot,var,iretcode,nfun)
2246 cdte       write(iout,*)'SUMSL MASK DIST return code is',iretcode,
2247 cdte     &                          ' eval ',nfun
2248 cdte       ieval=ieval+nfun
2249 cdte
2250 cdte       time1=MPI_WTIME()
2251 cdte       write (iout,'(a,f6.2,f8.2,a)')
2252 cdte     &        ' Time for mask dist min.',time1-time0,
2253 cdte     &         nfun/(time1-time0),'  eval/s'
2254 cdte       call flush(iout)
2255        if (debug) then
2256          call var_to_geom(nvar,var)
2257          call chainbuild
2258          call write_pdb(4000+in_pdb,'mask dist',etot)
2259        endif
2260 c
2261 c      switch off freezing of 2D and
2262 c      run full UNRES optimization with constrains
2263 c
2264        mask_r=.false.
2265 #ifdef MPI
2266        time0=MPI_WTIME()
2267 #else
2268        time0=tcpu()
2269 #endif
2270 cde       change=reduce(var)
2271 cde       if (check_var(var,info)) then
2272 cde         write(iout,*) 'error before dist'
2273 cde         call var_to_geom(nvar,var)
2274 cde         call chainbuild
2275 cde         call write_pdb(11000+in_pdb,'before dist',etot)
2276 cde       endif
2277
2278        call minimize(etot,var,iretcode,nfun)
2279
2280 cde        change=reduce(var)
2281 cde        if (check_var(var,info)) then
2282 cde          write(iout,*) 'error after dist',ico
2283 cde          call var_to_geom(nvar,var)
2284 cde          call chainbuild
2285 cde          call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
2286 cde        endif
2287        write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2288        ieval=ieval+nfun
2289
2290 #ifdef MPI
2291        time1=MPI_WTIME()
2292 #else
2293        time1=tcpu()
2294 #endif
2295        write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,
2296      &         nfun/(time1-time0),'  eval/s'
2297 cde       call etotal(energy(0))
2298 cde       write(iout,*) 'N7 after dist',energy(0)
2299        call flush(iout)
2300
2301        if (debug) then
2302         call var_to_geom(nvar,var)
2303         call chainbuild
2304         call write_pdb(in_pdb,linia,etot)
2305        endif
2306 c
2307 c      reset constrains
2308 c
2309        nhpb= nhpb0
2310        link_start=1
2311        link_end=nhpb
2312        wstrain=wstrain0
2313
2314        return
2315        end
2316 c--------------------------------------------------------
2317       subroutine softreg
2318       implicit real*8 (a-h,o-z)
2319       include 'DIMENSIONS'
2320 #ifdef MPI
2321       include 'mpif.h'
2322 #endif
2323       include 'COMMON.GEO'
2324       include 'COMMON.CHAIN'
2325       include 'COMMON.IOUNITS'
2326       include 'COMMON.VAR'
2327       include 'COMMON.CONTROL'
2328       include 'COMMON.SBRIDGE'
2329       include 'COMMON.FFIELD'
2330       include 'COMMON.MINIM'
2331       include 'COMMON.INTERACT'
2332 c
2333       include 'COMMON.DISTFIT'
2334       integer iff(maxres)
2335       double precision time0,time1
2336       double precision energy(0:n_ene),ee
2337       double precision var(maxvar)
2338       integer ieval
2339 c
2340       logical debug,ltest,fail
2341       character*50 linia
2342 c
2343       linia='test'
2344       debug=.true.
2345       in_pdb=0
2346
2347
2348
2349 c------------------------
2350 c
2351 c  freeze sec.elements
2352 c
2353        do i=1,nres
2354          mask_phi(i)=1
2355          mask_theta(i)=1
2356          mask_side(i)=1
2357          iff(i)=0
2358        enddo
2359
2360        do j=1,nbfrag
2361         do i=bfrag(1,j),bfrag(2,j)
2362          mask_phi(i)=0
2363          mask_theta(i)=0
2364          iff(i)=1
2365         enddo
2366         if (bfrag(3,j).le.bfrag(4,j)) then
2367          do i=bfrag(3,j),bfrag(4,j)
2368           mask_phi(i)=0
2369           mask_theta(i)=0
2370           iff(i)=1
2371          enddo
2372         else
2373          do i=bfrag(4,j),bfrag(3,j)
2374           mask_phi(i)=0
2375           mask_theta(i)=0
2376           iff(i)=1
2377          enddo
2378         endif
2379        enddo
2380        do j=1,nhfrag
2381         do i=hfrag(1,j),hfrag(2,j)
2382          mask_phi(i)=0
2383          mask_theta(i)=0
2384          iff(i)=1
2385         enddo
2386        enddo
2387        mask_r=.true.
2388
2389
2390
2391        nhpb0=nhpb
2392 c
2393 c store dist. constrains
2394 c
2395        do i=1,nres-3
2396          do j=i+3,nres
2397            if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2398             nhpb=nhpb+1
2399             ihpb(nhpb)=i
2400             jhpb(nhpb)=j
2401             forcon(nhpb)=0.1
2402             dhpb(nhpb)=DIST(i,j)
2403            endif
2404          enddo
2405        enddo
2406        call hpb_partition
2407
2408        if (debug) then
2409         call chainbuild
2410         call write_pdb(100+in_pdb,'input reg. structure',0d0)
2411        endif
2412
2413
2414        ipot0=ipot
2415        maxmin0=maxmin
2416        maxfun0=maxfun
2417        wstrain0=wstrain
2418        wang0=wang
2419 c
2420 c      run soft pot. optimization
2421 c
2422        ipot=6
2423        wang=3.0
2424        maxmin=2000
2425        maxfun=4000
2426        call geom_to_var(nvar,var)
2427 #ifdef MPI
2428        time0=MPI_WTIME()
2429 #else
2430        time0=tcpu()
2431 #endif
2432        call minimize(etot,var,iretcode,nfun)
2433
2434        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
2435 #ifdef MPI
2436        time1=MPI_WTIME()
2437 #else
2438        time1=tcpu()
2439 #endif
2440        write (iout,'(a,f6.2,f8.2,a)')'  Time for soft min.',time1-time0,
2441      &         nfun/(time1-time0),' SOFT eval/s'
2442        if (debug) then
2443          call var_to_geom(nvar,var)
2444          call chainbuild
2445          call write_pdb(300+in_pdb,'soft structure',etot)
2446        endif
2447 c
2448 c      run full UNRES optimization with constrains and frozen 2D
2449 c      the same variables as soft pot. optimizatio
2450 c
2451        ipot=ipot0
2452        wang=wang0
2453        maxmin=maxmin0
2454        maxfun=maxfun0
2455 #ifdef MPI
2456        time0=MPI_WTIME()
2457 #else
2458        time0=tcpu()
2459 #endif
2460        call minimize(etot,var,iretcode,nfun)
2461        write(iout,*)'SUMSL MASK DIST return code is',iretcode,
2462      &                          ' eval ',nfun
2463        ieval=nfun
2464
2465 #ifdef MPI
2466        time1=MPI_WTIME()
2467 #else
2468        time1=tcpu()
2469 #endif
2470        write (iout,'(a,f6.2,f8.2,a)')
2471      &        '  Time for mask dist min.',time1-time0,
2472      &         nfun/(time1-time0),'  eval/s'
2473        if (debug) then
2474          call var_to_geom(nvar,var)
2475          call chainbuild
2476          call write_pdb(400+in_pdb,'mask & dist',etot)
2477        endif
2478 c
2479 c      switch off constrains and
2480 c      run full UNRES optimization with frozen 2D
2481 c
2482
2483 c
2484 c      reset constrains
2485 c
2486        nhpb_c=nhpb
2487        nhpb=nhpb0
2488        link_start=1
2489        link_end=nhpb
2490        wstrain=wstrain0
2491
2492
2493 #ifdef MPI
2494        time0=MPI_WTIME()
2495 #else
2496        time0=tcpu()
2497 #endif
2498        call minimize(etot,var,iretcode,nfun)
2499        write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
2500        ieval=ieval+nfun
2501
2502 #ifdef MPI
2503        time1=MPI_WTIME()
2504 #else
2505        time1=tcpu()
2506 #endif
2507        write (iout,'(a,f6.2,f8.2,a)')'  Time for mask min.',time1-time0,
2508      &         nfun/(time1-time0),'  eval/s'
2509
2510
2511        if (debug) then
2512         call var_to_geom(nvar,var)
2513         call chainbuild
2514         call write_pdb(500+in_pdb,'mask 2d frozen',etot)
2515        endif
2516
2517        mask_r=.false.
2518
2519
2520 c
2521 c      run full UNRES optimization with constrains and NO frozen 2D
2522 c
2523
2524        nhpb=nhpb_c
2525        link_start=1
2526        link_end=nhpb
2527        maxfun=maxfun0/5
2528
2529        do ico=1,5
2530
2531        wstrain=wstrain0/ico
2532
2533 #ifdef MPI
2534        time0=MPI_WTIME()
2535 #else
2536        time0=tcpu()
2537 #endif
2538        call minimize(etot,var,iretcode,nfun)
2539        write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2540      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2541      &                          ' eval ',nfun
2542        ieval=nfun
2543
2544 #ifdef MPI
2545        time1=MPI_WTIME()
2546 #else
2547        time0=tcpu()
2548 #endif
2549        write (iout,'(a,f6.2,f8.2,a)')
2550      &        '  Time for dist min.',time1-time0,
2551      &         nfun/(time1-time0),'  eval/s'
2552        if (debug) then
2553          call var_to_geom(nvar,var)
2554          call chainbuild
2555          call write_pdb(600+in_pdb+ico,'dist cons',etot)
2556        endif
2557
2558        enddo
2559 c
2560        nhpb=nhpb0
2561        link_start=1
2562        link_end=nhpb
2563        wstrain=wstrain0
2564        maxfun=maxfun0
2565
2566
2567 c
2568       if (minim) then
2569 #ifdef MPI
2570        time0=MPI_WTIME()
2571 #else
2572        time0=tcpu()
2573 #endif
2574        call minimize(etot,var,iretcode,nfun)
2575        write(iout,*)'------------------------------------------------'
2576        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
2577      &  '+ DIST eval',ieval
2578
2579 #ifdef MPI
2580        time1=MPI_WTIME()
2581 #else
2582        time1=tcpu()
2583 #endif
2584        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
2585      &         nfun/(time1-time0),' eval/s'
2586
2587
2588        call var_to_geom(nvar,var)
2589        call chainbuild
2590        call write_pdb(999,'full min',etot)
2591       endif
2592
2593       return
2594       end
2595
2596
2597       subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
2598       implicit real*8 (a-h,o-z)
2599       include 'DIMENSIONS'
2600 #ifdef MPI
2601       include 'mpif.h'
2602 #endif
2603       include 'COMMON.GEO'
2604       include 'COMMON.VAR'
2605       include 'COMMON.INTERACT'
2606       include 'COMMON.IOUNITS'
2607       include 'COMMON.DISTFIT'
2608       include 'COMMON.SBRIDGE'
2609       include 'COMMON.CONTROL'
2610       include 'COMMON.FFIELD'
2611       include 'COMMON.MINIM'
2612       include 'COMMON.CHAIN'
2613       double precision time0,time1
2614       double precision energy(0:n_ene),ee
2615       double precision var(maxvar)
2616       integer jdata(5),isec(maxres)
2617 c
2618       jdata(1)=i1
2619       jdata(2)=i2
2620       jdata(3)=i3
2621       jdata(4)=i4
2622       jdata(5)=i5
2623
2624       call secondary2(.false.)
2625
2626       do i=1,nres
2627           isec(i)=0
2628       enddo
2629       do j=1,nbfrag
2630        do i=bfrag(1,j),bfrag(2,j)
2631           isec(i)=1
2632        enddo
2633        do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
2634           isec(i)=1
2635        enddo
2636       enddo
2637       do j=1,nhfrag
2638        do i=hfrag(1,j),hfrag(2,j)
2639           isec(i)=2
2640        enddo
2641       enddo
2642
2643 c
2644 c cut strands at the ends
2645 c
2646       if (jdata(2)-jdata(1).gt.3) then
2647        jdata(1)=jdata(1)+1
2648        jdata(2)=jdata(2)-1
2649        if (jdata(3).lt.jdata(4)) then
2650           jdata(3)=jdata(3)+1
2651           jdata(4)=jdata(4)-1
2652        else
2653           jdata(3)=jdata(3)-1
2654           jdata(4)=jdata(4)+1
2655        endif
2656       endif
2657
2658 cv      call chainbuild
2659 cv      call etotal(energy(0))
2660 cv      etot=energy(0)
2661 cv      write(iout,*) nnt,nct,etot
2662 cv      call write_pdb(ij*100,'first structure',etot)
2663 cv      write(iout,*) 'N16 test',(jdata(i),i=1,5)
2664
2665 c------------------------
2666 c      generate constrains
2667 c
2668        ishift=jdata(5)-2
2669        if(ishift.eq.0) ishift=-2
2670        nhpb0=nhpb
2671        call chainbuild
2672        do i=jdata(1),jdata(2)
2673         isec(i)=-1
2674         if(jdata(4).gt.jdata(3))then
2675          do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
2676             isec(j)=-1
2677 cd            print *,i,j,j+ishift
2678             nhpb=nhpb+1
2679             ihpb(nhpb)=i
2680             jhpb(nhpb)=j
2681             forcon(nhpb)=1000.0
2682             dhpb(nhpb)=DIST(i,j+ishift)
2683          enddo
2684         else
2685          do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
2686             isec(j)=-1
2687 cd            print *,i,j,j+ishift
2688             nhpb=nhpb+1
2689             ihpb(nhpb)=i
2690             jhpb(nhpb)=j
2691             forcon(nhpb)=1000.0
2692             dhpb(nhpb)=DIST(i,j+ishift)
2693          enddo
2694         endif
2695        enddo
2696
2697        do i=nnt,nct-2
2698          do j=i+2,nct
2699            if(isec(i).gt.0.or.isec(j).gt.0) then
2700 cd            print *,i,j
2701             nhpb=nhpb+1
2702             ihpb(nhpb)=i
2703             jhpb(nhpb)=j
2704             forcon(nhpb)=0.1
2705             dhpb(nhpb)=DIST(i,j)
2706            endif
2707          enddo
2708        enddo
2709
2710        call hpb_partition
2711
2712        call geom_to_var(nvar,var)
2713        maxfun0=maxfun
2714        wstrain0=wstrain
2715        maxfun=4000/5
2716
2717        do ico=1,5
2718
2719         wstrain=wstrain0/ico
2720
2721 cv        time0=MPI_WTIME()
2722         call minimize(etot,var,iretcode,nfun)
2723         write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2724      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2725      &                          ' eval ',nfun
2726         ieval=ieval+nfun
2727 cv        time1=MPI_WTIME()
2728 cv       write (iout,'(a,f6.2,f8.2,a)')
2729 cv     &        '  Time for dist min.',time1-time0,
2730 cv     &         nfun/(time1-time0),'  eval/s'
2731 cv         call var_to_geom(nvar,var)
2732 cv         call chainbuild
2733 cv         call write_pdb(ij*100+ico,'dist cons',etot)
2734
2735        enddo
2736 c
2737        nhpb=nhpb0
2738        call hpb_partition
2739        wstrain=wstrain0
2740        maxfun=maxfun0
2741 c
2742 cd      print *,etot
2743       wscloc0=wscloc
2744       wscloc=10.0
2745       call sc_move(nnt,nct,100,100d0,nft_sc,etot)
2746       wscloc=wscloc0
2747 cv      call chainbuild
2748 cv      call etotal(energy(0))
2749 cv      etot=energy(0)
2750 cv      call write_pdb(ij*100+10,'sc_move',etot)
2751 cd      call intout
2752 cd      print *,nft_sc,etot
2753
2754       return
2755       end
2756
2757       subroutine beta_zip(i1,i2,ieval,ij)
2758       implicit real*8 (a-h,o-z)
2759       include 'DIMENSIONS'
2760 #ifdef MPI
2761       include 'mpif.h'
2762 #endif
2763       include 'COMMON.GEO'
2764       include 'COMMON.VAR'
2765       include 'COMMON.INTERACT'
2766       include 'COMMON.IOUNITS'
2767       include 'COMMON.DISTFIT'
2768       include 'COMMON.SBRIDGE'
2769       include 'COMMON.CONTROL'
2770       include 'COMMON.FFIELD'
2771       include 'COMMON.MINIM'
2772       include 'COMMON.CHAIN'
2773       double precision time0,time1
2774       double precision energy(0:n_ene),ee
2775       double precision var(maxvar)
2776       character*10 test
2777
2778 cv      call chainbuild
2779 cv      call etotal(energy(0))
2780 cv      etot=energy(0)
2781 cv      write(test,'(2i5)') i1,i2
2782 cv      call write_pdb(ij*100,test,etot)
2783 cv      write(iout,*) 'N17 test',i1,i2,etot,ij
2784
2785 c
2786 c      generate constrains
2787 c
2788        nhpb0=nhpb
2789        nhpb=nhpb+1
2790        ihpb(nhpb)=i1
2791        jhpb(nhpb)=i2
2792        forcon(nhpb)=1000.0
2793        dhpb(nhpb)=4.0
2794
2795        call hpb_partition
2796
2797        call geom_to_var(nvar,var)
2798        maxfun0=maxfun
2799        wstrain0=wstrain
2800        maxfun=1000/5
2801
2802        do ico=1,5
2803         wstrain=wstrain0/ico
2804 cv        time0=MPI_WTIME()
2805         call minimize(etot,var,iretcode,nfun)
2806         write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2807      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2808      &                          ' eval ',nfun
2809         ieval=ieval+nfun
2810 cv        time1=MPI_WTIME()
2811 cv       write (iout,'(a,f6.2,f8.2,a)')
2812 cv     &        '  Time for dist min.',time1-time0,
2813 cv     &         nfun/(time1-time0),'  eval/s'
2814 c do not comment the next line
2815          call var_to_geom(nvar,var)
2816 cv         call chainbuild
2817 cv         call write_pdb(ij*100+ico,'dist cons',etot)
2818        enddo
2819
2820        nhpb=nhpb0
2821        call hpb_partition
2822        wstrain=wstrain0
2823        maxfun=maxfun0
2824
2825 cv      call etotal(energy(0))
2826 cv      etot=energy(0)
2827 cv      write(iout,*) 'N17 test end',i1,i2,etot,ij
2828
2829
2830       return
2831       end