source/unres/src_MD/src/test.F

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