source/unres/src_MD/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
 119
 120       subroutine test_local
 121       implicit real*8 (a-h,o-z)
 122       include 'DIMENSIONS'
 123       include 'COMMON.GEO'
 124       include 'COMMON.VAR'
 125       include 'COMMON.INTERACT'
 126       include 'COMMON.IOUNITS'
 127       double precision time0,time1
 128       double precision energy(0:n_ene),ee
 129       double precision varia(maxvar)
 130 c
 131       call chainbuild
 132 c      call geom_to_var(nvar,varia)
 133       call write_pdb(1,'first structure',0d0)
 134
 135       call etotal(energy(0))
 136       etot=energy(0)
 137       write(iout,*) nnt,nct,etot
 138
 139       write(iout,*) 'calling sc_move'
 140       call sc_move(nnt,nct,5,10d0,nft_sc,etot)
 141       write(iout,*) nft_sc,etot
 142       call write_pdb(2,'second structure',etot)
 143
 144       write(iout,*) 'calling local_move'
 145       call local_move_init(.false.)
 146       call local_move(24,29,20d0,50d0)
 147       call chainbuild
 148       call write_pdb(3,'third structure',etot)
 149
 150       write(iout,*) 'calling sc_move'
 151       call sc_move(24,29,5,10d0,nft_sc,etot)
 152       write(iout,*) nft_sc,etot
 153       call write_pdb(2,'last structure',etot)
 154
 155
 156       return
 157       end
 158
 159       subroutine test_sc
 160       implicit real*8 (a-h,o-z)
 161       include 'DIMENSIONS'
 162       include 'COMMON.GEO'
 163       include 'COMMON.VAR'
 164       include 'COMMON.INTERACT'
 165       include 'COMMON.IOUNITS'
 166       double precision time0,time1
 167       double precision energy(0:n_ene),ee
 168       double precision varia(maxvar)
 169 c
 170       call chainbuild
 171 c      call geom_to_var(nvar,varia)
 172       call write_pdb(1,'first structure',0d0)
 173
 174       call etotal(energy(0))
 175       etot=energy(0)
 176       write(iout,*) nnt,nct,etot
 177
 178       write(iout,*) 'calling sc_move'
 179
 180       call sc_move(nnt,nct,5,10d0,nft_sc,etot)
 181       write(iout,*) nft_sc,etot
 182       call write_pdb(2,'second structure',etot)
 183
 184       write(iout,*) 'calling sc_move 2nd time'
 185
 186       call sc_move(nnt,nct,5,1d0,nft_sc,etot)
 187       write(iout,*) nft_sc,etot
 188       call write_pdb(3,'last structure',etot)
 189       return
 190       end
 191 c--------------------------------------------------------
 192       subroutine bgrow(bstrand,nbstrand,in,ind,new)
 193       implicit real*8 (a-h,o-z)
 194       include 'DIMENSIONS'
 195       include 'COMMON.CHAIN'
 196       integer bstrand(maxres/3,6)
 197
 198       ishift=iabs(bstrand(in,ind+4)-new)
 199
 200       print *,'bgrow',bstrand(in,ind+4),new,ishift
 201
 202       bstrand(in,ind)=new
 203
 204       if(ind.eq.1)then
 205         bstrand(nbstrand,5)=bstrand(nbstrand,1)
 206         do i=1,nbstrand-1
 207           IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
 208           if (bstrand(i,5).lt.bstrand(i,6)) then
 209             bstrand(i,5)=bstrand(i,5)-ishift
 210           else
 211             bstrand(i,5)=bstrand(i,5)+ishift
 212           endif
 213           ENDIF
 214         enddo
 215       else
 216         bstrand(nbstrand,6)=bstrand(nbstrand,2)
 217         do i=1,nbstrand-1
 218           IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
 219           if (bstrand(i,6).lt.bstrand(i,5)) then
 220             bstrand(i,6)=bstrand(i,6)-ishift
 221           else
 222             bstrand(i,6)=bstrand(i,6)+ishift
 223           endif
 224           ENDIF
 225         enddo
 226       endif
 227
 228
 229       return
 230       end
 231
 232
 233 c-------------------------------------------------
 234
 235       subroutine secondary(lprint)
 236       implicit real*8 (a-h,o-z)
 237       include 'DIMENSIONS'
 238       include 'COMMON.CHAIN'
 239       include 'COMMON.IOUNITS'
 240       include 'COMMON.DISTFIT'
 241
 242       integer ncont,icont(2,maxres*maxres/2),isec(maxres,3)
 243       logical lprint,not_done
 244       real dcont(maxres*maxres/2),d
 245       real rcomp /7.0/
 246       real rbeta /5.2/
 247       real ralfa /5.2/
 248       real r310 /6.6/
 249       double precision xpi(3),xpj(3)
 250
 251
 252
 253       call chainbuild
 254 cd      call write_pdb(99,'sec structure',0d0)
 255       ncont=0
 256       nbfrag=0
 257       nhfrag=0
 258       do i=1,nres
 259         isec(i,1)=0
 260         isec(i,2)=0
 261         isec(i,3)=0
 262       enddo
 263
 264       do i=2,nres-3
 265         do k=1,3
 266           xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
 267         enddo
 268         do j=i+2,nres
 269           do k=1,3
 270              xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
 271           enddo
 272 cd        d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
 273 cd     &         (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
 274 cd     &         (c(3,i)-c(3,j))*(c(3,i)-c(3,j))
 275 cd          print *,'CA',i,j,d
 276           d =  (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) +
 277      &         (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) +
 278      &         (xpi(3)-xpj(3))*(xpi(3)-xpj(3))
 279           if ( d.lt.rcomp*rcomp) then
 280             ncont=ncont+1
 281             icont(1,ncont)=i
 282             icont(2,ncont)=j
 283             dcont(ncont)=sqrt(d)
 284           endif
 285         enddo
 286       enddo
 287       if (lprint) then
 288         write (iout,*)
 289         write (iout,'(a)') '#PP contact map distances:'
 290         do i=1,ncont
 291           write (iout,'(3i4,f10.5)')
 292      &     i,icont(1,i),icont(2,i),dcont(i)
 293         enddo
 294       endif
 295
 296 c finding parallel beta
 297 cd      write (iout,*) '------- looking for parallel beta -----------'
 298       nbeta=0
 299       nstrand=0
 300       do i=1,ncont
 301         i1=icont(1,i)
 302         j1=icont(2,i)
 303         if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and.
 304      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
 305      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and.
 306      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and.
 307      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and.
 308      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
 309      &     ) then
 310           ii1=i1
 311           jj1=j1
 312 cd         write (iout,*) i1,j1,dcont(i)
 313           not_done=.true.
 314           do while (not_done)
 315            i1=i1+1
 316            j1=j1+1
 317             do j=1,ncont
 318               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)
 319      &              .and. dcont(j).le.rbeta .and.
 320      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
 321      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and.
 322      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and.
 323      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and.
 324      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
 325      &                            ) goto 5
 326             enddo
 327             not_done=.false.
 328   5         continue
 329 cd            write (iout,*) i1,j1,dcont(j),not_done
 330           enddo
 331           j1=j1-1
 332           i1=i1-1
 333           if (i1-ii1.gt.1) then
 334             ii1=max0(ii1-1,1)
 335             jj1=max0(jj1-1,1)
 336             nbeta=nbeta+1
 337             if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
 338
 339             nbfrag=nbfrag+1
 340             bfrag(1,nbfrag)=ii1
 341             bfrag(2,nbfrag)=i1
 342             bfrag(3,nbfrag)=jj1
 343             bfrag(4,nbfrag)=j1
 344
 345             do ij=ii1,i1
 346              isec(ij,1)=isec(ij,1)+1
 347              isec(ij,1+isec(ij,1))=nbeta
 348             enddo
 349             do ij=jj1,j1
 350              isec(ij,1)=isec(ij,1)+1
 351              isec(ij,1+isec(ij,1))=nbeta
 352             enddo
 353
 354            if(lprint) then
 355             nstrand=nstrand+1
 356             if (nbeta.le.9) then
 357               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 358      &          "DefPropRes 'strand",nstrand,
 359      &          "' 'num = ",ii1-1,"..",i1-1,"'"
 360             else
 361               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 362      &          "DefPropRes 'strand",nstrand,
 363      &          "' 'num = ",ii1-1,"..",i1-1,"'"
 364             endif
 365             nstrand=nstrand+1
 366             if (nbeta.le.9) then
 367               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 368      &          "DefPropRes 'strand",nstrand,
 369      &          "' 'num = ",jj1-1,"..",j1-1,"'"
 370             else
 371               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 372      &          "DefPropRes 'strand",nstrand,
 373      &          "' 'num = ",jj1-1,"..",j1-1,"'"
 374             endif
 375               write(12,'(a8,4i4)')
 376      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
 377            endif
 378           endif
 379         endif
 380       enddo
 381
 382 c finding antiparallel beta
 383 cd      write (iout,*) '--------- looking for antiparallel beta ---------'
 384
 385       do i=1,ncont
 386         i1=icont(1,i)
 387         j1=icont(2,i)
 388         if (dcont(i).le.rbeta.and.
 389      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
 390      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and.
 391      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and.
 392      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and.
 393      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
 394      &     ) then
 395           ii1=i1
 396           jj1=j1
 397 cd          write (iout,*) i1,j1,dcont(i)
 398
 399           not_done=.true.
 400           do while (not_done)
 401            i1=i1+1
 402            j1=j1-1
 403             do j=1,ncont
 404               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
 405      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
 406      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and.
 407      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and.
 408      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and.
 409      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
 410      &           .and. dcont(j).le.rbeta ) goto 6
 411             enddo
 412             not_done=.false.
 413   6         continue
 414 cd            write (iout,*) i1,j1,dcont(j),not_done
 415           enddo
 416           i1=i1-1
 417           j1=j1+1
 418           if (i1-ii1.gt.1) then
 419             if(lprint)write (iout,*)'antiparallel beta',
 420      &                   nbeta,ii1-1,i1,jj1,j1-1
 421
 422             nbfrag=nbfrag+1
 423             bfrag(1,nbfrag)=max0(ii1-1,1)
 424             bfrag(2,nbfrag)=i1
 425             bfrag(3,nbfrag)=jj1
 426             bfrag(4,nbfrag)=max0(j1-1,1)
 427
 428             nbeta=nbeta+1
 429             iii1=max0(ii1-1,1)
 430             do ij=iii1,i1
 431              isec(ij,1)=isec(ij,1)+1
 432              isec(ij,1+isec(ij,1))=nbeta
 433             enddo
 434             jjj1=max0(j1-1,1)
 435             do ij=jjj1,jj1
 436              isec(ij,1)=isec(ij,1)+1
 437              isec(ij,1+isec(ij,1))=nbeta
 438             enddo
 439
 440
 441            if (lprint) then
 442             nstrand=nstrand+1
 443             if (nstrand.le.9) then
 444               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 445      &          "DefPropRes 'strand",nstrand,
 446      &          "' 'num = ",ii1-2,"..",i1-1,"'"
 447             else
 448               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 449      &          "DefPropRes 'strand",nstrand,
 450      &          "' 'num = ",ii1-2,"..",i1-1,"'"
 451             endif
 452             nstrand=nstrand+1
 453             if (nstrand.le.9) then
 454               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 455      &          "DefPropRes 'strand",nstrand,
 456      &          "' 'num = ",j1-2,"..",jj1-1,"'"
 457             else
 458               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 459      &          "DefPropRes 'strand",nstrand,
 460      &          "' 'num = ",j1-2,"..",jj1-1,"'"
 461             endif
 462               write(12,'(a8,4i4)')
 463      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
 464            endif
 465           endif
 466         endif
 467       enddo
 468
 469       if (nstrand.gt.0.and.lprint) then
 470         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
 471         do i=2,nstrand
 472          if (i.le.9) then
 473           write(12,'(a9,i1,$)') " | strand",i
 474          else
 475           write(12,'(a9,i2,$)') " | strand",i
 476          endif
 477         enddo
 478         write(12,'(a1)') "'"
 479       endif
 480
 481
 482 c finding alpha or 310 helix
 483
 484       nhelix=0
 485       do i=1,ncont
 486         i1=icont(1,i)
 487         j1=icont(2,i)
 488         if (j1.eq.i1+3.and.dcont(i).le.r310
 489      &     .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
 490 cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
 491 cd          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
 492           ii1=i1
 493           jj1=j1
 494           if (isec(ii1,1).eq.0) then
 495             not_done=.true.
 496           else
 497             not_done=.false.
 498           endif
 499           do while (not_done)
 500             i1=i1+1
 501             j1=j1+1
 502             do j=1,ncont
 503               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
 504             enddo
 505             not_done=.false.
 506   10        continue
 507 cd            write (iout,*) i1,j1,not_done
 508           enddo
 509           j1=j1-1
 510           if (j1-ii1.gt.4) then
 511             nhelix=nhelix+1
 512 cd            write (iout,*)'helix',nhelix,ii1,j1
 513
 514             nhfrag=nhfrag+1
 515             hfrag(1,nhfrag)=ii1
 516             hfrag(2,nhfrag)=max0(j1-1,1)
 517
 518             do ij=ii1,j1
 519              isec(ij,1)=-1
 520             enddo
 521            if (lprint) then
 522             write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
 523             if (nhelix.le.9) then
 524               write(12,'(a17,i1,a9,i3,a2,i3,a1)')
 525      &          "DefPropRes 'helix",nhelix,
 526      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 527             else
 528               write(12,'(a17,i2,a9,i3,a2,i3,a1)')
 529      &          "DefPropRes 'helix",nhelix,
 530      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 531             endif
 532            endif
 533           endif
 534         endif
 535       enddo
 536
 537       if (nhelix.gt.0.and.lprint) then
 538         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
 539         do i=2,nhelix
 540          if (nhelix.le.9) then
 541           write(12,'(a8,i1,$)') " | helix",i
 542          else
 543           write(12,'(a8,i2,$)') " | helix",i
 544          endif
 545         enddo
 546         write(12,'(a1)') "'"
 547       endif
 548
 549       if (lprint) then
 550        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
 551        write(12,'(a20)') "XMacStand ribbon.mac"
 552       endif
 553
 554
 555       return
 556       end
 557 c----------------------------------------------------------------------------
 558
 559       subroutine write_pdb(npdb,titelloc,ee)
 560       implicit real*8 (a-h,o-z)
 561       include 'DIMENSIONS'
 562       include 'COMMON.IOUNITS'
 563       character*50 titelloc1
 564       character*(*) titelloc
 565       character*3 zahl
 566       character*5 liczba5
 567       double precision ee
 568       integer npdb,ilen
 569       external ilen
 570
 571       titelloc1=titelloc
 572       lenpre=ilen(prefix)
 573       if (npdb.lt.1000) then
 574        call numstr(npdb,zahl)
 575        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
 576       else
 577         if (npdb.lt.10000) then
 578          write(liczba5,'(i1,i4)') 0,npdb
 579         else
 580          write(liczba5,'(i5)') npdb
 581         endif
 582         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
 583       endif
 584       call pdbout(ee,titelloc1,ipdb)
 585       close(ipdb)
 586       return
 587       end
 588
 589 c--------------------------------------------------------
 590       subroutine softreg
 591       implicit real*8 (a-h,o-z)
 592       include 'DIMENSIONS'
 593 #ifdef MPI
 594       include 'mpif.h'
 595 #endif
 596       include 'COMMON.GEO'
 597       include 'COMMON.CHAIN'
 598       include 'COMMON.IOUNITS'
 599       include 'COMMON.VAR'
 600       include 'COMMON.CONTROL'
 601       include 'COMMON.SBRIDGE'
 602       include 'COMMON.FFIELD'
 603       include 'COMMON.MINIM'
 604       include 'COMMON.INTERACT'
 605 c
 606       include 'COMMON.DISTFIT'
 607       integer iff(maxres)
 608       double precision time0,time1
 609       double precision energy(0:n_ene),ee
 610       double precision var(maxvar)
 611       integer ieval
 612 c
 613       logical debug,ltest,fail
 614       character*50 linia
 615 c
 616       linia='test'
 617       debug=.true.
 618       in_pdb=0
 619
 620
 621
 622 c------------------------
 623 c
 624 c  freeze sec.elements
 625 c
 626        do i=1,nres
 627          mask_phi(i)=1
 628          mask_theta(i)=1
 629          mask_side(i)=1
 630          iff(i)=0
 631        enddo
 632
 633        do j=1,nbfrag
 634         do i=bfrag(1,j),bfrag(2,j)
 635          mask_phi(i)=0
 636          mask_theta(i)=0
 637          iff(i)=1
 638         enddo
 639         if (bfrag(3,j).le.bfrag(4,j)) then
 640          do i=bfrag(3,j),bfrag(4,j)
 641           mask_phi(i)=0
 642           mask_theta(i)=0
 643           iff(i)=1
 644          enddo
 645         else
 646          do i=bfrag(4,j),bfrag(3,j)
 647           mask_phi(i)=0
 648           mask_theta(i)=0
 649           iff(i)=1
 650          enddo
 651         endif
 652        enddo
 653        do j=1,nhfrag
 654         do i=hfrag(1,j),hfrag(2,j)
 655          mask_phi(i)=0
 656          mask_theta(i)=0
 657          iff(i)=1
 658         enddo
 659        enddo
 660        mask_r=.true.
 661
 662
 663
 664        nhpb0=nhpb
 665 c
 666 c store dist. constrains
 667 c
 668        do i=1,nres-3
 669          do j=i+3,nres
 670            if ( iff(i).eq.1.and.iff(j).eq.1 ) then
 671             nhpb=nhpb+1
 672             ihpb(nhpb)=i
 673             jhpb(nhpb)=j
 674             forcon(nhpb)=0.1
 675             dhpb(nhpb)=DIST(i,j)
 676            endif
 677          enddo
 678        enddo
 679        call hpb_partition
 680
 681        if (debug) then
 682         call chainbuild
 683         call write_pdb(100+in_pdb,'input reg. structure',0d0)
 684        endif
 685
 686
 687        ipot0=ipot
 688        maxmin0=maxmin
 689        maxfun0=maxfun
 690        wstrain0=wstrain
 691        wang0=wang
 692 c
 693 c      run soft pot. optimization
 694 c
 695        ipot=6
 696        wang=3.0
 697        maxmin=2000
 698        maxfun=4000
 699        call geom_to_var(nvar,var)
 700 #ifdef MPI
 701        time0=MPI_WTIME()
 702 #else
 703        time0=tcpu()
 704 #endif
 705        call minimize(etot,var,iretcode,nfun)
 706
 707        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun
 708 #ifdef MPI
 709        time1=MPI_WTIME()
 710 #else
 711        time1=tcpu()
 712 #endif
 713        write (iout,'(a,f6.2,f8.2,a)')'  Time for soft min.',time1-time0,
 714      &         nfun/(time1-time0),' SOFT eval/s'
 715        if (debug) then
 716          call var_to_geom(nvar,var)
 717          call chainbuild
 718          call write_pdb(300+in_pdb,'soft structure',etot)
 719        endif
 720 c
 721 c      run full UNRES optimization with constrains and frozen 2D
 722 c      the same variables as soft pot. optimizatio
 723 c
 724        ipot=ipot0
 725        wang=wang0
 726        maxmin=maxmin0
 727        maxfun=maxfun0
 728 #ifdef MPI
 729        time0=MPI_WTIME()
 730 #else
 731        time0=tcpu()
 732 #endif
 733        call minimize(etot,var,iretcode,nfun)
 734        write(iout,*)'SUMSL MASK DIST return code is',iretcode,
 735      &                          ' eval ',nfun
 736        ieval=nfun
 737 #ifdef MPI
 738        time1=MPI_WTIME()
 739 #else
 740        time1=tcpu()
 741 #endif
 742        write (iout,'(a,f6.2,f8.2,a)')
 743      &        '  Time for mask dist min.',time1-time0,
 744      &         nfun/(time1-time0),'  eval/s'
 745        if (debug) then
 746          call var_to_geom(nvar,var)
 747          call chainbuild
 748          call write_pdb(400+in_pdb,'mask & dist',etot)
 749        endif
 750 c
 751 c      switch off constrains and
 752 c      run full UNRES optimization with frozen 2D
 753 c
 754
 755 c
 756 c      reset constrains
 757 c
 758        nhpb_c=nhpb
 759        nhpb=nhpb0
 760        link_start=1
 761        link_end=nhpb
 762        wstrain=wstrain0
 763
 764 #ifdef MPI
 765        time0=MPI_WTIME()
 766 #else
 767        time0=tcpu()
 768 #endif
 769        call minimize(etot,var,iretcode,nfun)
 770        write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
 771        ieval=ieval+nfun
 772 #ifdef MPI
 773        time1=MPI_WTIME()
 774 #else
 775        time1=tcpu()
 776 #endif
 777        write (iout,'(a,f6.2,f8.2,a)')'  Time for mask min.',time1-time0,
 778      &         nfun/(time1-time0),'  eval/s'
 779
 780
 781        if (debug) then
 782         call var_to_geom(nvar,var)
 783         call chainbuild
 784         call write_pdb(500+in_pdb,'mask 2d frozen',etot)
 785        endif
 786
 787        mask_r=.false.
 788
 789
 790 c
 791 c      run full UNRES optimization with constrains and NO frozen 2D
 792 c
 793
 794        nhpb=nhpb_c
 795        link_start=1
 796        link_end=nhpb
 797        maxfun=maxfun0/5
 798
 799        do ico=1,5
 800
 801        wstrain=wstrain0/ico
 802 #ifdef MPI
 803        time0=MPI_WTIME()
 804 #else
 805        time0=tcpu()
 806 #endif
 807        call minimize(etot,var,iretcode,nfun)
 808        write(iout,'(a10,f6.3,a14,i3,a6,i5)')
 809      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
 810      &                          ' eval ',nfun
 811        ieval=nfun
 812 #ifdef MPI
 813        time1=MPI_WTIME()
 814 #else
 815        time1=tcpu()
 816 #endif
 817        write (iout,'(a,f6.2,f8.2,a)')
 818      &        '  Time for dist min.',time1-time0,
 819      &         nfun/(time1-time0),'  eval/s'
 820        if (debug) then
 821          call var_to_geom(nvar,var)
 822          call chainbuild
 823          call write_pdb(600+in_pdb+ico,'dist cons',etot)
 824        endif
 825
 826        enddo
 827 c
 828        nhpb=nhpb0
 829        link_start=1
 830        link_end=nhpb
 831        wstrain=wstrain0
 832        maxfun=maxfun0
 833
 834
 835 c
 836       if (minim) then
 837 #ifdef MPI
 838        time0=MPI_WTIME()
 839 #else
 840        time0=tcpu()
 841 #endif
 842        call minimize(etot,var,iretcode,nfun)
 843        write(iout,*)'------------------------------------------------'
 844        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
 845      &  '+ DIST eval',ieval
 846 #ifdef MPI
 847        time1=MPI_WTIME()
 848 #else
 849        time1=tcpu()
 850 #endif
 851        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
 852      &         nfun/(time1-time0),' eval/s'
 853
 854
 855        call var_to_geom(nvar,var)
 856        call chainbuild
 857        call write_pdb(999,'full min',etot)
 858       endif
 859
 860       return
 861       end
 862
 863