source/unres/src-HCD-5D/contact_cp.F

   1       subroutine contact_cp(var,var2,iff,ieval,in_pdb)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.SBRIDGE'
   5       include 'COMMON.FFIELD'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.FRAG'
   8       include 'COMMON.VAR'
   9       include 'COMMON.CHAIN'
  10       include 'COMMON.MINIM'
  11
  12       character*50 linia
  13       integer nf,ij(4)
  14       double precision energy(0:n_ene)
  15       double precision var(maxvar),var2(maxvar)
  16       double precision time0,time1
  17       integer iff(maxres),ieval
  18       double precision theta1(maxres),phi1(maxres),alph1(maxres),
  19      &                 omeg1(maxres)
  20       logical debug
  21
  22       debug=.false.
  23 c      debug=.true.
  24       if (ieval.eq.-1) debug=.true.
  25
  26
  27 c
  28 c store selected dist. constrains from 1st structure
  29 c
  30 #ifdef OSF
  31 c     Intercept NaNs in the coordinates
  32 c      write(iout,*) (var(i),i=1,nvar)
  33       x_sum=0.D0
  34       do i=1,nvar
  35         x_sum=x_sum+var(i)
  36       enddo
  37       if (x_sum.ne.x_sum) then
  38         write(iout,*)" *** contact_cp : Found NaN in coordinates"
  39         call flush(iout)
  40         print *," *** contact_cp : Found NaN in coordinates"
  41         return
  42       endif
  43 #endif
  44
  45
  46        call var_to_geom(nvar,var)
  47        call chainbuild
  48        nhpb0=nhpb
  49        ind=0
  50        do i=1,nres-3
  51          do j=i+3,nres
  52           ind=ind+1
  53           if ( iff(i).eq.1.and.iff(j).eq.1 ) then
  54 c            d0(ind)=DIST(i,j)
  55 c            w(ind)=10.0
  56             nhpb=nhpb+1
  57             ihpb(nhpb)=i
  58             jhpb(nhpb)=j
  59             forcon(nhpb)=10.0
  60             dhpb(nhpb)=DIST(i,j)
  61           else
  62 c            w(ind)=0.0
  63           endif
  64          enddo
  65        enddo
  66        call hpb_partition
  67
  68        do i=1,nres
  69         theta1(i)=theta(i)
  70         phi1(i)=phi(i)
  71         alph1(i)=alph(i)
  72         omeg1(i)=omeg(i)
  73        enddo
  74
  75 c
  76 c  freeze sec.elements from 2nd structure
  77 c
  78        do i=1,nres
  79          mask_phi(i)=1
  80          mask_theta(i)=1
  81          mask_side(i)=1
  82        enddo
  83
  84        call var_to_geom(nvar,var2)
  85        call secondary2(debug)
  86        do j=1,nbfrag
  87         do i=bfrag(1,j),bfrag(2,j)
  88 c         mask(i)=0
  89          mask_phi(i)=0
  90          mask_theta(i)=0
  91         enddo
  92         if (bfrag(3,j).le.bfrag(4,j)) then
  93          do i=bfrag(3,j),bfrag(4,j)
  94 c          mask(i)=0
  95           mask_phi(i)=0
  96           mask_theta(i)=0
  97          enddo
  98         else
  99          do i=bfrag(4,j),bfrag(3,j)
 100 c          mask(i)=0
 101           mask_phi(i)=0
 102           mask_theta(i)=0
 103          enddo
 104         endif
 105        enddo
 106        do j=1,nhfrag
 107         do i=hfrag(1,j),hfrag(2,j)
 108 c         mask(i)=0
 109          mask_phi(i)=0
 110          mask_theta(i)=0
 111         enddo
 112        enddo
 113        mask_r=.true.
 114
 115 c
 116 c      copy selected res from 1st to 2nd structure
 117 c
 118
 119        do i=1,nres
 120           if ( iff(i).eq.1 ) then
 121               theta(i)=theta1(i)
 122               phi(i)=phi1(i)
 123               alph(i)=alph1(i)
 124               omeg(i)=omeg1(i)
 125           endif
 126        enddo
 127
 128       if(debug) then
 129 c
 130 c     prepare description in linia variable
 131 c
 132         iwsk=0
 133         nf=0
 134         if (iff(1).eq.1) then
 135           iwsk=1
 136           nf=nf+1
 137           ij(nf)=1
 138         endif
 139         do i=2,nres
 140            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
 141              iwsk=1
 142              nf=nf+1
 143              ij(nf)=i
 144            endif
 145            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
 146              iwsk=0
 147              nf=nf+1
 148              ij(nf)=i-1
 149            endif
 150         enddo
 151         if (iff(nres).eq.1) then
 152           nf=nf+1
 153           ij(nf)=nres
 154         endif
 155
 156         write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
 157      &                     "SELECT",ij(1)-1,"-",ij(2)-1,
 158      &                         ",",ij(3)-1,"-",ij(4)-1
 159
 160       endif
 161 c
 162 c     run optimization
 163 c
 164       call contact_cp_min(var,ieval,in_pdb,linia,debug)
 165
 166       return
 167       end
 168
 169       subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
 170 c
 171 c    input : theta,phi,alph,omeg,in_pdb,linia,debug
 172 c    output : var,ieval
 173 c
 174       implicit real*8 (a-h,o-z)
 175       include 'DIMENSIONS'
 176 #ifdef MPI
 177       include 'mpif.h'
 178 #endif
 179       include 'COMMON.SBRIDGE'
 180       include 'COMMON.FFIELD'
 181       include 'COMMON.IOUNITS'
 182       include 'COMMON.FRAG'
 183       include 'COMMON.VAR'
 184       include 'COMMON.CHAIN'
 185       include 'COMMON.MINIM'
 186
 187       character*50 linia
 188       integer nf,ij(4)
 189       double precision energy(0:n_ene)
 190       double precision var(maxvar)
 191       double precision time0,time1
 192       integer ieval,info(3)
 193       logical debug,fail,check_var,reduce,change
 194
 195        write(iout,'(a20,i6,a20)')
 196      &             '------------------',in_pdb,'-------------------'
 197
 198        if (debug) then
 199         call chainbuild
 200         call write_pdb(1000+in_pdb,'combined structure',0d0)
 201 #ifdef MPI
 202         time0=MPI_WTIME()
 203 #else
 204         time0=tcpu()
 205 #endif
 206        endif
 207
 208 c
 209 c     run optimization of distances
 210 c
 211 c     uses d0(),w() and mask() for frozen 2D
 212 c
 213 ctest---------------------------------------------
 214 ctest       NX=NRES-3
 215 ctest       NY=((NRES-4)*(NRES-5))/2
 216 ctest       call distfit(debug,5000)
 217
 218        do i=1,nres
 219          mask_side(i)=0
 220        enddo
 221
 222        ipot01=ipot
 223        maxmin01=maxmin
 224        maxfun01=maxfun
 225 c       wstrain01=wstrain
 226        wsc01=wsc
 227        wscp01=wscp
 228        welec01=welec
 229        wvdwpp01=wvdwpp
 230 c      wang01=wang
 231        wscloc01=wscloc
 232        wtor01=wtor
 233        wtor_d01=wtor_d
 234
 235        ipot=6
 236        maxmin=2000
 237        maxfun=4000
 238 c       wstrain=1.0
 239        wsc=0.0
 240        wscp=0.0
 241        welec=0.0
 242        wvdwpp=0.0
 243 c      wang=0.0
 244        wscloc=0.0
 245        wtor=0.0
 246        wtor_d=0.0
 247
 248        call geom_to_var(nvar,var)
 249 cde       change=reduce(var)
 250        if (check_var(var,info)) then
 251           write(iout,*) 'cp_min error in input'
 252           print *,'cp_min error in input'
 253           return
 254        endif
 255
 256 cd       call etotal(energy(0))
 257 cd       call enerprint(energy(0))
 258 cd       call check_eint
 259
 260 #ifdef MPI
 261        time0=MPI_WTIME()
 262 #else
 263        time0=tcpu()
 264 #endif
 265 cdtest       call minimize(etot,var,iretcode,nfun)
 266 cdtest       write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun
 267 #ifdef MPI
 268        time1=MPI_WTIME()
 269 #else
 270        time1=tcpu()
 271 #endif
 272
 273 cd       call etotal(energy(0))
 274 cd       call enerprint(energy(0))
 275 cd       call check_eint
 276
 277        do i=1,nres
 278          mask_side(i)=1
 279        enddo
 280
 281        ipot=ipot01
 282        maxmin=maxmin01
 283        maxfun=maxfun01
 284 c       wstrain=wstrain01
 285        wsc=wsc01
 286        wscp=wscp01
 287        welec=welec01
 288        wvdwpp=wvdwpp01
 289 c      wang=wang01
 290        wscloc=wscloc01
 291        wtor=wtor01
 292        wtor_d=wtor_d01
 293 ctest--------------------------------------------------
 294
 295        if(debug) then
 296 #ifdef MPI
 297         time1=MPI_WTIME()
 298 #else
 299         time1=tcpu()
 300 #endif
 301         write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
 302         call write_pdb(2000+in_pdb,'distfit structure',0d0)
 303        endif
 304
 305
 306        ipot0=ipot
 307        maxmin0=maxmin
 308        maxfun0=maxfun
 309        wstrain0=wstrain
 310 c
 311 c      run soft pot. optimization
 312 c         with constrains:
 313 c             nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition
 314 c         and frozen 2D:
 315 c             mask_phi(),mask_theta(),mask_side(),mask_r
 316 c
 317        ipot=6
 318        maxmin=2000
 319        maxfun=4000
 320
 321 cde       change=reduce(var)
 322 cde       if (check_var(var,info)) write(iout,*) 'error before soft'
 323 #ifdef MPI
 324        time0=MPI_WTIME()
 325 #else
 326        time0=tcpu()
 327 #endif
 328        call minimize(etot,var,iretcode,nfun)
 329
 330        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
 331 #ifdef MPI
 332        time1=MPI_WTIME()
 333 #else
 334        time1=tcpu()
 335 #endif
 336        write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
 337      &         nfun/(time1-time0),' SOFT eval/s'
 338        if (debug) then
 339          call var_to_geom(nvar,var)
 340          call chainbuild
 341          call write_pdb(3000+in_pdb,'soft structure',etot)
 342        endif
 343 c
 344 c      run full UNRES optimization with constrains and frozen 2D
 345 c      the same variables as soft pot. optimizatio
 346 c
 347        ipot=ipot0
 348        maxmin=maxmin0
 349        maxfun=maxfun0
 350 c
 351 c check overlaps before calling full UNRES minim
 352 c
 353        call var_to_geom(nvar,var)
 354        call chainbuild
 355        call etotal(energy(0))
 356 #ifdef OSF
 357        write(iout,*) 'N7 ',energy(0)
 358        if (energy(0).ne.energy(0)) then
 359         write(iout,*) 'N7 error - gives NaN',energy(0)
 360        endif
 361 #endif
 362        ieval=1
 363        if (energy(1).eq.1.0d20) then
 364          write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
 365          call overlap_sc(fail)
 366          if(.not.fail) then
 367            call etotal(energy(0))
 368            ieval=ieval+1
 369            write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
 370          else
 371            mask_r=.false.
 372            nhpb= nhpb0
 373            link_start=1
 374            link_end=nhpb
 375            wstrain=wstrain0
 376            return
 377          endif
 378        endif
 379        call flush(iout)
 380 c
 381 cdte       time0=MPI_WTIME()
 382 cde       change=reduce(var)
 383 cde       if (check_var(var,info)) then
 384 cde         write(iout,*) 'error before mask dist'
 385 cde         call var_to_geom(nvar,var)
 386 cde         call chainbuild
 387 cde         call write_pdb(10000+in_pdb,'before mask dist',etot)
 388 cde       endif
 389 cdte       call minimize(etot,var,iretcode,nfun)
 390 cdte       write(iout,*)'SUMSL MASK DIST return code is',iretcode,
 391 cdte     &                          ' eval ',nfun
 392 cdte       ieval=ieval+nfun
 393 cdte
 394 cdte       time1=MPI_WTIME()
 395 cdte       write (iout,'(a,f6.2,f8.2,a)')
 396 cdte     &        ' Time for mask dist min.',time1-time0,
 397 cdte     &         nfun/(time1-time0),'  eval/s'
 398 cdte       call flush(iout)
 399        if (debug) then
 400          call var_to_geom(nvar,var)
 401          call chainbuild
 402          call write_pdb(4000+in_pdb,'mask dist',etot)
 403        endif
 404 c
 405 c      switch off freezing of 2D and
 406 c      run full UNRES optimization with constrains
 407 c
 408        mask_r=.false.
 409 #ifdef MPI
 410        time0=MPI_WTIME()
 411 #else
 412        time0=tcpu()
 413 #endif
 414 cde       change=reduce(var)
 415 cde       if (check_var(var,info)) then
 416 cde         write(iout,*) 'error before dist'
 417 cde         call var_to_geom(nvar,var)
 418 cde         call chainbuild
 419 cde         call write_pdb(11000+in_pdb,'before dist',etot)
 420 cde       endif
 421
 422        call minimize(etot,var,iretcode,nfun)
 423
 424 cde        change=reduce(var)
 425 cde        if (check_var(var,info)) then
 426 cde          write(iout,*) 'error after dist',ico
 427 cde          call var_to_geom(nvar,var)
 428 cde          call chainbuild
 429 cde          call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
 430 cde        endif
 431        write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
 432        ieval=ieval+nfun
 433
 434 #ifdef MPI
 435        time1=MPI_WTIME()
 436 #else
 437        time1=tcpu()
 438 #endif
 439        write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,
 440      &         nfun/(time1-time0),'  eval/s'
 441 cde       call etotal(energy(0))
 442 cde       write(iout,*) 'N7 after dist',energy(0)
 443 c       call flush(iout)
 444
 445        if (debug) then
 446         call var_to_geom(nvar,var)
 447         call chainbuild
 448         call write_pdb(in_pdb,linia,etot)
 449        endif
 450 c
 451 c      reset constrains
 452 c
 453        nhpb= nhpb0
 454        link_start=1
 455        link_end=nhpb
 456        wstrain=wstrain0
 457
 458        return
 459        end
 460 c--------------------------------------------------------
 461       subroutine softreg
 462       implicit real*8 (a-h,o-z)
 463       include 'DIMENSIONS'
 464 #ifdef MPI
 465       include 'mpif.h'
 466 #endif
 467       include 'COMMON.GEO'
 468       include 'COMMON.CHAIN'
 469       include 'COMMON.IOUNITS'
 470       include 'COMMON.VAR'
 471       include 'COMMON.CONTROL'
 472       include 'COMMON.SBRIDGE'
 473       include 'COMMON.FFIELD'
 474       include 'COMMON.MINIM'
 475       include 'COMMON.INTERACT'
 476 c
 477       include 'COMMON.FRAG'
 478       integer iff(maxres)
 479       double precision time0,time1
 480       double precision energy(0:n_ene),ee
 481       double precision var(maxvar)
 482       integer ieval
 483 c
 484       logical debug,ltest,fail
 485       character*50 linia
 486 c
 487       linia='test'
 488       debug=.true.
 489       in_pdb=0
 490
 491
 492
 493 c------------------------
 494 c
 495 c  freeze sec.elements
 496 c
 497        do i=1,nres
 498          mask_phi(i)=1
 499          mask_theta(i)=1
 500          mask_side(i)=1
 501          iff(i)=0
 502        enddo
 503
 504        do j=1,nbfrag
 505         do i=bfrag(1,j),bfrag(2,j)
 506          mask_phi(i)=0
 507          mask_theta(i)=0
 508          iff(i)=1
 509         enddo
 510         if (bfrag(3,j).le.bfrag(4,j)) then
 511          do i=bfrag(3,j),bfrag(4,j)
 512           mask_phi(i)=0
 513           mask_theta(i)=0
 514           iff(i)=1
 515          enddo
 516         else
 517          do i=bfrag(4,j),bfrag(3,j)
 518           mask_phi(i)=0
 519           mask_theta(i)=0
 520           iff(i)=1
 521          enddo
 522         endif
 523        enddo
 524        do j=1,nhfrag
 525         do i=hfrag(1,j),hfrag(2,j)
 526          mask_phi(i)=0
 527          mask_theta(i)=0
 528          iff(i)=1
 529         enddo
 530        enddo
 531        mask_r=.true.
 532
 533
 534
 535        nhpb0=nhpb
 536 c
 537 c store dist. constrains
 538 c
 539        do i=1,nres-3
 540          do j=i+3,nres
 541            if ( iff(i).eq.1.and.iff(j).eq.1 ) then
 542             nhpb=nhpb+1
 543             ihpb(nhpb)=i
 544             jhpb(nhpb)=j
 545             forcon(nhpb)=0.1
 546             dhpb(nhpb)=DIST(i,j)
 547            endif
 548          enddo
 549        enddo
 550        call hpb_partition
 551
 552        if (debug) then
 553         call chainbuild
 554         call write_pdb(100+in_pdb,'input reg. structure',0d0)
 555        endif
 556
 557
 558        ipot0=ipot
 559        maxmin0=maxmin
 560        maxfun0=maxfun
 561        wstrain0=wstrain
 562        wang0=wang
 563 c
 564 c      run soft pot. optimization
 565 c
 566        ipot=6
 567        wang=3.0
 568        maxmin=2000
 569        maxfun=4000
 570        call geom_to_var(nvar,var)
 571 #ifdef MPI
 572        time0=MPI_WTIME()
 573 #else
 574        time0=tcpu()
 575 #endif
 576        call minimize(etot,var,iretcode,nfun)
 577
 578        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
 579 #ifdef MPI
 580        time1=MPI_WTIME()
 581 #else
 582        time1=tcpu()
 583 #endif
 584        write (iout,'(a,f6.2,f8.2,a)')'  Time for soft min.',time1-time0,
 585      &         nfun/(time1-time0),' SOFT eval/s'
 586        if (debug) then
 587          call var_to_geom(nvar,var)
 588          call chainbuild
 589          call write_pdb(300+in_pdb,'soft structure',etot)
 590        endif
 591 c
 592 c      run full UNRES optimization with constrains and frozen 2D
 593 c      the same variables as soft pot. optimizatio
 594 c
 595        ipot=ipot0
 596        wang=wang0
 597        maxmin=maxmin0
 598        maxfun=maxfun0
 599 #ifdef MPI
 600        time0=MPI_WTIME()
 601 #else
 602        time0=tcpu()
 603 #endif
 604        call minimize(etot,var,iretcode,nfun)
 605        write(iout,*)'SUMSL MASK DIST return code is',iretcode,
 606      &                          ' eval ',nfun
 607        ieval=nfun
 608
 609 #ifdef MPI
 610        time1=MPI_WTIME()
 611 #else
 612        time1=tcpu()
 613 #endif
 614        write (iout,'(a,f6.2,f8.2,a)')
 615      &        '  Time for mask dist min.',time1-time0,
 616      &         nfun/(time1-time0),'  eval/s'
 617        if (debug) then
 618          call var_to_geom(nvar,var)
 619          call chainbuild
 620          call write_pdb(400+in_pdb,'mask & dist',etot)
 621        endif
 622 c
 623 c      switch off constrains and
 624 c      run full UNRES optimization with frozen 2D
 625 c
 626
 627 c
 628 c      reset constrains
 629 c
 630        nhpb_c=nhpb
 631        nhpb=nhpb0
 632        link_start=1
 633        link_end=nhpb
 634        wstrain=wstrain0
 635
 636
 637 #ifdef MPI
 638        time0=MPI_WTIME()
 639 #else
 640        time0=tcpu()
 641 #endif
 642        call minimize(etot,var,iretcode,nfun)
 643        write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
 644        ieval=ieval+nfun
 645
 646 #ifdef MPI
 647        time1=MPI_WTIME()
 648 #else
 649        time1=tcpu()
 650 #endif
 651        write (iout,'(a,f6.2,f8.2,a)')'  Time for mask min.',time1-time0,
 652      &         nfun/(time1-time0),'  eval/s'
 653
 654
 655        if (debug) then
 656         call var_to_geom(nvar,var)
 657         call chainbuild
 658         call write_pdb(500+in_pdb,'mask 2d frozen',etot)
 659        endif
 660
 661        mask_r=.false.
 662
 663
 664 c
 665 c      run full UNRES optimization with constrains and NO frozen 2D
 666 c
 667
 668        nhpb=nhpb_c
 669        link_start=1
 670        link_end=nhpb
 671        maxfun=maxfun0/5
 672
 673        do ico=1,5
 674
 675        wstrain=wstrain0/ico
 676
 677 #ifdef MPI
 678        time0=MPI_WTIME()
 679 #else
 680        time0=tcpu()
 681 #endif
 682        call minimize(etot,var,iretcode,nfun)
 683        write(iout,'(a10,f6.3,a14,i3,a6,i5)')
 684      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
 685      &                          ' eval ',nfun
 686        ieval=nfun
 687
 688 #ifdef MPI
 689        time1=MPI_WTIME()
 690 #else
 691        time0=tcpu()
 692 #endif
 693        write (iout,'(a,f6.2,f8.2,a)')
 694      &        '  Time for dist min.',time1-time0,
 695      &         nfun/(time1-time0),'  eval/s'
 696        if (debug) then
 697          call var_to_geom(nvar,var)
 698          call chainbuild
 699          call write_pdb(600+in_pdb+ico,'dist cons',etot)
 700        endif
 701
 702        enddo
 703 c
 704        nhpb=nhpb0
 705        link_start=1
 706        link_end=nhpb
 707        wstrain=wstrain0
 708        maxfun=maxfun0
 709
 710
 711 c
 712       if (minim) then
 713 #ifdef MPI
 714        time0=MPI_WTIME()
 715 #else
 716        time0=tcpu()
 717 #endif
 718        call minimize(etot,var,iretcode,nfun)
 719        write(iout,*)'------------------------------------------------'
 720        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
 721      &  '+ DIST eval',ieval
 722
 723 #ifdef MPI
 724        time1=MPI_WTIME()
 725 #else
 726        time1=tcpu()
 727 #endif
 728        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
 729      &         nfun/(time1-time0),' eval/s'
 730
 731
 732        call var_to_geom(nvar,var)
 733        call chainbuild
 734        call write_pdb(999,'full min',etot)
 735       endif
 736
 737       return
 738       end
 739
 740
 741       subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
 742       implicit real*8 (a-h,o-z)
 743       include 'DIMENSIONS'
 744 #ifdef MPI
 745       include 'mpif.h'
 746 #endif
 747       include 'COMMON.GEO'
 748       include 'COMMON.VAR'
 749       include 'COMMON.INTERACT'
 750       include 'COMMON.IOUNITS'
 751       include 'COMMON.FRAG'
 752       include 'COMMON.SBRIDGE'
 753       include 'COMMON.CONTROL'
 754       include 'COMMON.FFIELD'
 755       include 'COMMON.MINIM'
 756       include 'COMMON.CHAIN'
 757       double precision time0,time1
 758       double precision energy(0:n_ene),ee
 759       double precision var(maxvar)
 760       integer jdata(5),isec(maxres)
 761 c
 762       jdata(1)=i1
 763       jdata(2)=i2
 764       jdata(3)=i3
 765       jdata(4)=i4
 766       jdata(5)=i5
 767
 768       call secondary2(.false.)
 769
 770       do i=1,nres
 771           isec(i)=0
 772       enddo
 773       do j=1,nbfrag
 774        do i=bfrag(1,j),bfrag(2,j)
 775           isec(i)=1
 776        enddo
 777        do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
 778           isec(i)=1
 779        enddo
 780       enddo
 781       do j=1,nhfrag
 782        do i=hfrag(1,j),hfrag(2,j)
 783           isec(i)=2
 784        enddo
 785       enddo
 786
 787 c
 788 c cut strands at the ends
 789 c
 790       if (jdata(2)-jdata(1).gt.3) then
 791        jdata(1)=jdata(1)+1
 792        jdata(2)=jdata(2)-1
 793        if (jdata(3).lt.jdata(4)) then
 794           jdata(3)=jdata(3)+1
 795           jdata(4)=jdata(4)-1
 796        else
 797           jdata(3)=jdata(3)-1
 798           jdata(4)=jdata(4)+1
 799        endif
 800       endif
 801
 802 cv      call chainbuild
 803 cv      call etotal(energy(0))
 804 cv      etot=energy(0)
 805 cv      write(iout,*) nnt,nct,etot
 806 cv      call write_pdb(ij*100,'first structure',etot)
 807 cv      write(iout,*) 'N16 test',(jdata(i),i=1,5)
 808
 809 c------------------------
 810 c      generate constrains
 811 c
 812        ishift=jdata(5)-2
 813        if(ishift.eq.0) ishift=-2
 814        nhpb0=nhpb
 815        call chainbuild
 816        do i=jdata(1),jdata(2)
 817         isec(i)=-1
 818         if(jdata(4).gt.jdata(3))then
 819          do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
 820             isec(j)=-1
 821 cd            print *,i,j,j+ishift
 822             nhpb=nhpb+1
 823             ihpb(nhpb)=i
 824             jhpb(nhpb)=j
 825             forcon(nhpb)=1000.0
 826             dhpb(nhpb)=DIST(i,j+ishift)
 827          enddo
 828         else
 829          do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
 830             isec(j)=-1
 831 cd            print *,i,j,j+ishift
 832             nhpb=nhpb+1
 833             ihpb(nhpb)=i
 834             jhpb(nhpb)=j
 835             forcon(nhpb)=1000.0
 836             dhpb(nhpb)=DIST(i,j+ishift)
 837          enddo
 838         endif
 839        enddo
 840
 841        do i=nnt,nct-2
 842          do j=i+2,nct
 843            if(isec(i).gt.0.or.isec(j).gt.0) then
 844 cd            print *,i,j
 845             nhpb=nhpb+1
 846             ihpb(nhpb)=i
 847             jhpb(nhpb)=j
 848             forcon(nhpb)=0.1
 849             dhpb(nhpb)=DIST(i,j)
 850            endif
 851          enddo
 852        enddo
 853
 854        call hpb_partition
 855
 856        call geom_to_var(nvar,var)
 857        maxfun0=maxfun
 858        wstrain0=wstrain
 859        maxfun=4000/5
 860
 861        do ico=1,5
 862
 863         wstrain=wstrain0/ico
 864
 865 cv        time0=MPI_WTIME()
 866         call minimize(etot,var,iretcode,nfun)
 867         write(iout,'(a10,f6.3,a14,i3,a6,i5)')
 868      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
 869      &                          ' eval ',nfun
 870         ieval=ieval+nfun
 871 cv        time1=MPI_WTIME()
 872 cv       write (iout,'(a,f6.2,f8.2,a)')
 873 cv     &        '  Time for dist min.',time1-time0,
 874 cv     &         nfun/(time1-time0),'  eval/s'
 875 cv         call var_to_geom(nvar,var)
 876 cv         call chainbuild
 877 cv         call write_pdb(ij*100+ico,'dist cons',etot)
 878
 879        enddo
 880 c
 881        nhpb=nhpb0
 882        call hpb_partition
 883        wstrain=wstrain0
 884        maxfun=maxfun0
 885 c
 886 cd      print *,etot
 887       wscloc0=wscloc
 888       wscloc=10.0
 889       call sc_move(nnt,nct,100,100d0,nft_sc,etot)
 890       wscloc=wscloc0
 891 cv      call chainbuild
 892 cv      call etotal(energy(0))
 893 cv      etot=energy(0)
 894 cv      call write_pdb(ij*100+10,'sc_move',etot)
 895 cd      call intout
 896 cd      print *,nft_sc,etot
 897
 898       return
 899       end
 900
 901       subroutine beta_zip(i1,i2,ieval,ij)
 902       implicit real*8 (a-h,o-z)
 903       include 'DIMENSIONS'
 904 #ifdef MPI
 905       include 'mpif.h'
 906 #endif
 907       include 'COMMON.GEO'
 908       include 'COMMON.VAR'
 909       include 'COMMON.INTERACT'
 910       include 'COMMON.IOUNITS'
 911       include 'COMMON.FRAG'
 912       include 'COMMON.SBRIDGE'
 913       include 'COMMON.CONTROL'
 914       include 'COMMON.FFIELD'
 915       include 'COMMON.MINIM'
 916       include 'COMMON.CHAIN'
 917       double precision time0,time1
 918       double precision energy(0:n_ene),ee
 919       double precision var(maxvar)
 920       character*10 test
 921
 922 cv      call chainbuild
 923 cv      call etotal(energy(0))
 924 cv      etot=energy(0)
 925 cv      write(test,'(2i5)') i1,i2
 926 cv      call write_pdb(ij*100,test,etot)
 927 cv      write(iout,*) 'N17 test',i1,i2,etot,ij
 928
 929 c
 930 c      generate constrains
 931 c
 932        nhpb0=nhpb
 933        nhpb=nhpb+1
 934        ihpb(nhpb)=i1
 935        jhpb(nhpb)=i2
 936        forcon(nhpb)=1000.0
 937        dhpb(nhpb)=4.0
 938
 939        call hpb_partition
 940
 941        call geom_to_var(nvar,var)
 942        maxfun0=maxfun
 943        wstrain0=wstrain
 944        maxfun=1000/5
 945
 946        do ico=1,5
 947         wstrain=wstrain0/ico
 948 cv        time0=MPI_WTIME()
 949         call minimize(etot,var,iretcode,nfun)
 950         write(iout,'(a10,f6.3,a14,i3,a6,i5)')
 951      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
 952      &                          ' eval ',nfun
 953         ieval=ieval+nfun
 954 cv        time1=MPI_WTIME()
 955 cv       write (iout,'(a,f6.2,f8.2,a)')
 956 cv     &        '  Time for dist min.',time1-time0,
 957 cv     &         nfun/(time1-time0),'  eval/s'
 958 c do not comment the next line
 959          call var_to_geom(nvar,var)
 960 cv         call chainbuild
 961 cv         call write_pdb(ij*100+ico,'dist cons',etot)
 962        enddo
 963
 964        nhpb=nhpb0
 965        call hpb_partition
 966        wstrain=wstrain0
 967        maxfun=maxfun0
 968
 969 cv      call etotal(energy(0))
 970 cv      etot=energy(0)
 971 cv      write(iout,*) 'N17 test end',i1,i2,etot,ij
 972
 973
 974       return
 975       end
 976 c----------------------------------------------------------------------------
 977
 978       subroutine write_pdb(npdb,titelloc,ee)
 979       implicit real*8 (a-h,o-z)
 980       include 'DIMENSIONS'
 981       include 'COMMON.IOUNITS'
 982       character*50 titelloc1
 983       character*(*) titelloc
 984       character*3 zahl
 985       character*5 liczba5
 986       double precision ee
 987       integer npdb,ilen
 988       external ilen
 989
 990       titelloc1=titelloc
 991       lenpre=ilen(prefix)
 992       if (npdb.lt.1000) then
 993        call numstr(npdb,zahl)
 994        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
 995       else
 996         if (npdb.lt.10000) then
 997          write(liczba5,'(i1,i4)') 0,npdb
 998         else
 999          write(liczba5,'(i5)') npdb
1000         endif
1001         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
1002       endif
1003       call pdbout(ee,titelloc1,ipdb)
1004       close(ipdb)
1005       return
1006       end
1007