source/unres/src-HCD-5D/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 c      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 c      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 c      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 c      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 c      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 c      include 'COMMON.DISTFIT'
1363
1364       integer ncont,icont(2,maxres*maxint_res),isec(maxres,3)
1365       logical lprint,not_done
1366       real dcont(maxres*maxint_res),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-----------------------------------------------------------