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