Backup before merge with adasko
[unres.git] / source / unres / src_MD_DFA / test.F
diff --git a/source/unres/src_MD_DFA/test.F b/source/unres/src_MD_DFA/test.F
deleted file mode 100644 (file)
index 0140ee5..0000000
+++ /dev/null
@@ -1,863 +0,0 @@
-      subroutine test
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-#ifdef MPI
-      include 'mpif.h'
-#endif
-      include 'COMMON.GEO'
-      include 'COMMON.VAR'
-      include 'COMMON.INTERACT'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.DISTFIT'
-      include 'COMMON.SBRIDGE'
-      include 'COMMON.CONTROL'
-      include 'COMMON.FFIELD'
-      include 'COMMON.MINIM'
-      include 'COMMON.CHAIN'
-      double precision time0,time1
-      double precision energy(0:n_ene),ee
-      double precision var(maxvar),var1(maxvar)
-      integer j1,j2
-      logical debug,accepted
-      debug=.true.
-      
-
-      call geom_to_var(nvar,var1)
-      call chainbuild
-      call etotal(energy(0))
-      etot=energy(0)
-      call rmsd(rms)
-      write(iout,*) 'etot=',0,etot,rms
-      call secondary2(.false.)
-
-      call write_pdb(0,'first structure',etot)
-
-      j1=13
-      j2=21
-      da=180.0*deg2rad
-
-
-
-       temp=3000.0d0
-       betbol=1.0D0/(1.9858D-3*temp)
-       jr=iran_num(j1,j2)
-       d=ran_number(-pi,pi)
-c       phi(jr)=pinorm(phi(jr)+d)
-       call chainbuild
-       call etotal(energy(0))
-       etot0=energy(0)
-       call rmsd(rms)
-       write(iout,*) 'etot=',1,etot0,rms
-       call write_pdb(1,'perturb structure',etot0)
-
-      do i=2,500,2
-       jr=iran_num(j1,j2)
-       d=ran_number(-da,da)
-       phiold=phi(jr)
-       phi(jr)=pinorm(phi(jr)+d)
-       call chainbuild
-       call etotal(energy(0))
-       etot=energy(0)
-
-       if (etot.lt.etot0) then 
-          accepted=.true.
-       else
-          accepted=.false.
-          xxr=ran_number(0.0D0,1.0D0)
-          xxh=betbol*(etot-etot0)
-          if (xxh.lt.50.0D0) then
-            xxh=dexp(-xxh)
-            if (xxh.gt.xxr) accepted=.true. 
-          endif
-       endif
-       accepted=.true.
-c       print *,etot0,etot,accepted
-       if (accepted) then 
-          etot0=etot
-          call rmsd(rms)
-          write(iout,*) 'etot=',i,etot,rms
-          call write_pdb(i,'MC structure',etot)
-c minimize
-c        call geom_to_var(nvar,var1)
-        call sc_move(2,nres-1,1,10d0,nft_sc,etot)
-        call geom_to_var(nvar,var)
-        call minimize(etot,var,iretcode,nfun)
-        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
-        call var_to_geom(nvar,var)
-        call chainbuild
-        call rmsd(rms)
-        write(iout,*) 'etot mcm=',i,etot,rms
-        call write_pdb(i+1,'MCM structure',etot)
-        call var_to_geom(nvar,var1)
-c --------
-       else
-          phi(jr)=phiold
-       endif
-      enddo
-
-c minimize
-c       call sc_move(2,nres-1,1,10d0,nft_sc,etot)
-c       call geom_to_var(nvar,var)
-c
-c       call chainbuild        
-c       call write_pdb(998 ,'sc min',etot)
-c
-c       call minimize(etot,var,iretcode,nfun)
-c       write(iout,*)'------------------------------------------------'
-c       write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
-c      
-c       call var_to_geom(nvar,var)
-c       call chainbuild        
-c       call write_pdb(999,'full min',etot)
-
-
-      return
-      end
-
-
-
-
-      subroutine test_local
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'COMMON.GEO'
-      include 'COMMON.VAR'
-      include 'COMMON.INTERACT'
-      include 'COMMON.IOUNITS'
-      double precision time0,time1
-      double precision energy(0:n_ene),ee
-      double precision varia(maxvar)
-c
-      call chainbuild
-c      call geom_to_var(nvar,varia)
-      call write_pdb(1,'first structure',0d0)
-
-      call etotal(energy(0))
-      etot=energy(0)
-      write(iout,*) nnt,nct,etot
-
-      write(iout,*) 'calling sc_move'
-      call sc_move(nnt,nct,5,10d0,nft_sc,etot)
-      write(iout,*) nft_sc,etot
-      call write_pdb(2,'second structure',etot)
-
-      write(iout,*) 'calling local_move'
-      call local_move_init(.false.)
-      call local_move(24,29,20d0,50d0)     
-      call chainbuild
-      call write_pdb(3,'third structure',etot)
-
-      write(iout,*) 'calling sc_move'
-      call sc_move(24,29,5,10d0,nft_sc,etot)
-      write(iout,*) nft_sc,etot
-      call write_pdb(2,'last structure',etot)
-
-
-      return
-      end
-
-      subroutine test_sc
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'COMMON.GEO'
-      include 'COMMON.VAR'
-      include 'COMMON.INTERACT'
-      include 'COMMON.IOUNITS'
-      double precision time0,time1
-      double precision energy(0:n_ene),ee
-      double precision varia(maxvar)
-c
-      call chainbuild
-c      call geom_to_var(nvar,varia)
-      call write_pdb(1,'first structure',0d0)
-
-      call etotal(energy(0))
-      etot=energy(0)
-      write(iout,*) nnt,nct,etot
-
-      write(iout,*) 'calling sc_move'
-
-      call sc_move(nnt,nct,5,10d0,nft_sc,etot)
-      write(iout,*) nft_sc,etot
-      call write_pdb(2,'second structure',etot)
-
-      write(iout,*) 'calling sc_move 2nd time'
-
-      call sc_move(nnt,nct,5,1d0,nft_sc,etot)
-      write(iout,*) nft_sc,etot
-      call write_pdb(3,'last structure',etot)
-      return
-      end
-c--------------------------------------------------------
-      subroutine bgrow(bstrand,nbstrand,in,ind,new)
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'COMMON.CHAIN'
-      integer bstrand(maxres/3,6)
-
-      ishift=iabs(bstrand(in,ind+4)-new)
-
-      print *,'bgrow',bstrand(in,ind+4),new,ishift
-
-      bstrand(in,ind)=new
-
-      if(ind.eq.1)then
-        bstrand(nbstrand,5)=bstrand(nbstrand,1)
-        do i=1,nbstrand-1
-          IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
-          if (bstrand(i,5).lt.bstrand(i,6)) then 
-            bstrand(i,5)=bstrand(i,5)-ishift
-          else
-            bstrand(i,5)=bstrand(i,5)+ishift
-          endif
-          ENDIF
-        enddo
-      else
-        bstrand(nbstrand,6)=bstrand(nbstrand,2)
-        do i=1,nbstrand-1
-          IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
-          if (bstrand(i,6).lt.bstrand(i,5)) then 
-            bstrand(i,6)=bstrand(i,6)-ishift
-          else
-            bstrand(i,6)=bstrand(i,6)+ishift
-          endif
-          ENDIF
-        enddo
-      endif
-
-
-      return
-      end
-
-
-c-------------------------------------------------
-
-      subroutine secondary(lprint)
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'COMMON.CHAIN'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.DISTFIT'
-
-      integer ncont,icont(2,maxres*maxres/2),isec(maxres,3)
-      logical lprint,not_done
-      real dcont(maxres*maxres/2),d
-      real rcomp /7.0/ 
-      real rbeta /5.2/
-      real ralfa /5.2/
-      real r310 /6.6/
-      double precision xpi(3),xpj(3)
-
-
-
-      call chainbuild
-cd      call write_pdb(99,'sec structure',0d0)
-      ncont=0
-      nbfrag=0
-      nhfrag=0
-      do i=1,nres
-        isec(i,1)=0
-        isec(i,2)=0
-        isec(i,3)=0
-      enddo
-
-      do i=2,nres-3
-        do k=1,3        
-          xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
-        enddo
-        do j=i+2,nres
-          do k=1,3
-             xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
-          enddo
-cd       d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
-cd     &         (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
-cd     &         (c(3,i)-c(3,j))*(c(3,i)-c(3,j)) 
-cd          print *,'CA',i,j,d
-          d =  (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) +
-     &         (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) +
-     &         (xpi(3)-xpj(3))*(xpi(3)-xpj(3)) 
-         if ( d.lt.rcomp*rcomp) then
-            ncont=ncont+1
-            icont(1,ncont)=i
-            icont(2,ncont)=j
-            dcont(ncont)=sqrt(d)
-          endif
-        enddo
-      enddo
-      if (lprint) then
-        write (iout,*)
-        write (iout,'(a)') '#PP contact map distances:'
-        do i=1,ncont
-          write (iout,'(3i4,f10.5)') 
-     &     i,icont(1,i),icont(2,i),dcont(i) 
-        enddo
-      endif
-
-c finding parallel beta
-cd      write (iout,*) '------- looking for parallel beta -----------'
-      nbeta=0
-      nstrand=0
-      do i=1,ncont
-        i1=icont(1,i)
-        j1=icont(2,i)
-        if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and.
-     &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
-     &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. 
-     &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. 
-     &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. 
-     &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
-     &     ) then
-          ii1=i1
-          jj1=j1
-cd         write (iout,*) i1,j1,dcont(i)
-          not_done=.true.
-          do while (not_done)
-           i1=i1+1
-           j1=j1+1
-            do j=1,ncont
-              if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)
-     &              .and. dcont(j).le.rbeta .and.
-     &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
-     &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. 
-     &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. 
-     &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. 
-     &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
-     &                            ) goto 5
-            enddo
-            not_done=.false.
-  5         continue
-cd            write (iout,*) i1,j1,dcont(j),not_done
-          enddo
-          j1=j1-1
-          i1=i1-1
-          if (i1-ii1.gt.1) then
-            ii1=max0(ii1-1,1)
-            jj1=max0(jj1-1,1)
-            nbeta=nbeta+1
-            if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
-
-            nbfrag=nbfrag+1
-            bfrag(1,nbfrag)=ii1
-            bfrag(2,nbfrag)=i1
-            bfrag(3,nbfrag)=jj1
-            bfrag(4,nbfrag)=j1 
-
-            do ij=ii1,i1
-             isec(ij,1)=isec(ij,1)+1
-             isec(ij,1+isec(ij,1))=nbeta
-            enddo
-            do ij=jj1,j1
-             isec(ij,1)=isec(ij,1)+1
-             isec(ij,1+isec(ij,1))=nbeta
-            enddo
-
-           if(lprint) then 
-            nstrand=nstrand+1
-            if (nbeta.le.9) then
-              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'strand",nstrand,
-     &          "' 'num = ",ii1-1,"..",i1-1,"'"
-            else
-              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'strand",nstrand,
-     &          "' 'num = ",ii1-1,"..",i1-1,"'"
-            endif
-            nstrand=nstrand+1
-            if (nbeta.le.9) then
-              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'strand",nstrand,
-     &          "' 'num = ",jj1-1,"..",j1-1,"'"
-            else
-              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'strand",nstrand,
-     &          "' 'num = ",jj1-1,"..",j1-1,"'"
-            endif
-              write(12,'(a8,4i4)')
-     &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
-           endif
-          endif
-        endif
-      enddo
-
-c finding antiparallel beta
-cd      write (iout,*) '--------- looking for antiparallel beta ---------'
-
-      do i=1,ncont
-        i1=icont(1,i)
-        j1=icont(2,i)
-        if (dcont(i).le.rbeta.and.
-     &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
-     &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. 
-     &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. 
-     &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. 
-     &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
-     &     ) then
-          ii1=i1
-          jj1=j1
-cd          write (iout,*) i1,j1,dcont(i)
-
-          not_done=.true.
-          do while (not_done)
-           i1=i1+1
-           j1=j1-1
-            do j=1,ncont
-              if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
-     &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
-     &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. 
-     &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. 
-     &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. 
-     &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
-     &           .and. dcont(j).le.rbeta ) goto 6
-            enddo
-            not_done=.false.
-  6         continue
-cd            write (iout,*) i1,j1,dcont(j),not_done
-          enddo
-          i1=i1-1
-          j1=j1+1
-          if (i1-ii1.gt.1) then
-            if(lprint)write (iout,*)'antiparallel beta',
-     &                   nbeta,ii1-1,i1,jj1,j1-1
-
-            nbfrag=nbfrag+1
-            bfrag(1,nbfrag)=max0(ii1-1,1)
-            bfrag(2,nbfrag)=i1
-            bfrag(3,nbfrag)=jj1
-            bfrag(4,nbfrag)=max0(j1-1,1) 
-
-            nbeta=nbeta+1
-            iii1=max0(ii1-1,1)
-            do ij=iii1,i1
-             isec(ij,1)=isec(ij,1)+1
-             isec(ij,1+isec(ij,1))=nbeta
-            enddo
-            jjj1=max0(j1-1,1)  
-            do ij=jjj1,jj1
-             isec(ij,1)=isec(ij,1)+1
-             isec(ij,1+isec(ij,1))=nbeta
-            enddo
-
-
-           if (lprint) then
-            nstrand=nstrand+1
-            if (nstrand.le.9) then
-              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'strand",nstrand,
-     &          "' 'num = ",ii1-2,"..",i1-1,"'"
-            else
-              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'strand",nstrand,
-     &          "' 'num = ",ii1-2,"..",i1-1,"'"
-            endif
-            nstrand=nstrand+1
-            if (nstrand.le.9) then
-              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'strand",nstrand,
-     &          "' 'num = ",j1-2,"..",jj1-1,"'"
-            else
-              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'strand",nstrand,
-     &          "' 'num = ",j1-2,"..",jj1-1,"'"
-            endif
-              write(12,'(a8,4i4)')
-     &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
-           endif
-          endif
-        endif
-      enddo
-
-      if (nstrand.gt.0.and.lprint) then
-        write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
-        do i=2,nstrand
-         if (i.le.9) then
-          write(12,'(a9,i1,$)') " | strand",i
-         else
-          write(12,'(a9,i2,$)') " | strand",i
-         endif
-        enddo
-        write(12,'(a1)') "'"
-      endif
-
-       
-c finding alpha or 310 helix
-
-      nhelix=0
-      do i=1,ncont
-        i1=icont(1,i)
-        j1=icont(2,i)
-        if (j1.eq.i1+3.and.dcont(i).le.r310
-     &     .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
-cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
-cd          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
-          ii1=i1
-          jj1=j1
-          if (isec(ii1,1).eq.0) then 
-            not_done=.true.
-          else
-            not_done=.false.
-          endif
-          do while (not_done)
-            i1=i1+1
-            j1=j1+1
-            do j=1,ncont
-              if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
-            enddo
-            not_done=.false.
-  10        continue
-cd            write (iout,*) i1,j1,not_done
-          enddo
-          j1=j1-1
-          if (j1-ii1.gt.4) then
-            nhelix=nhelix+1
-cd            write (iout,*)'helix',nhelix,ii1,j1
-
-            nhfrag=nhfrag+1
-            hfrag(1,nhfrag)=ii1
-            hfrag(2,nhfrag)=max0(j1-1,1)
-
-            do ij=ii1,j1
-             isec(ij,1)=-1
-            enddo
-           if (lprint) then
-            write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
-            if (nhelix.le.9) then
-              write(12,'(a17,i1,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'helix",nhelix,
-     &          "' 'num = ",ii1-1,"..",j1-2,"'"
-            else
-              write(12,'(a17,i2,a9,i3,a2,i3,a1)') 
-     &          "DefPropRes 'helix",nhelix,
-     &          "' 'num = ",ii1-1,"..",j1-2,"'"
-            endif
-           endif
-          endif
-        endif
-      enddo
-       
-      if (nhelix.gt.0.and.lprint) then
-        write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
-        do i=2,nhelix
-         if (nhelix.le.9) then
-          write(12,'(a8,i1,$)') " | helix",i
-         else
-          write(12,'(a8,i2,$)') " | helix",i
-         endif
-        enddo
-        write(12,'(a1)') "'"
-      endif
-
-      if (lprint) then
-       write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
-       write(12,'(a20)') "XMacStand ribbon.mac"
-      endif
-
-
-      return
-      end
-c----------------------------------------------------------------------------
-
-      subroutine write_pdb(npdb,titelloc,ee)
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'COMMON.IOUNITS'
-      character*50 titelloc1                                                     
-      character*(*) titelloc
-      character*3 zahl   
-      character*5 liczba5
-      double precision ee
-      integer npdb,ilen
-      external ilen
-
-      titelloc1=titelloc
-      lenpre=ilen(prefix)
-      if (npdb.lt.1000) then
-       call numstr(npdb,zahl)
-       open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
-      else
-        if (npdb.lt.10000) then                              
-         write(liczba5,'(i1,i4)') 0,npdb
-        else   
-         write(liczba5,'(i5)') npdb
-        endif
-        open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
-      endif
-      call pdbout(ee,titelloc1,ipdb)
-      close(ipdb)
-      return
-      end
-
-c--------------------------------------------------------
-      subroutine softreg
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-#ifdef MPI
-      include 'mpif.h'
-#endif
-      include 'COMMON.GEO'
-      include 'COMMON.CHAIN'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.VAR'
-      include 'COMMON.CONTROL'
-      include 'COMMON.SBRIDGE'
-      include 'COMMON.FFIELD'
-      include 'COMMON.MINIM'
-      include 'COMMON.INTERACT'
-c
-      include 'COMMON.DISTFIT'       
-      integer iff(maxres)
-      double precision time0,time1
-      double precision energy(0:n_ene),ee
-      double precision var(maxvar)
-      integer ieval
-c
-      logical debug,ltest,fail
-      character*50 linia
-c
-      linia='test'
-      debug=.true.
-      in_pdb=0
-
-
-
-c------------------------
-c
-c  freeze sec.elements 
-c
-       do i=1,nres
-         mask_phi(i)=1
-         mask_theta(i)=1
-         mask_side(i)=1
-         iff(i)=0
-       enddo
-
-       do j=1,nbfrag
-        do i=bfrag(1,j),bfrag(2,j)
-         mask_phi(i)=0
-         mask_theta(i)=0
-         iff(i)=1
-        enddo
-        if (bfrag(3,j).le.bfrag(4,j)) then 
-         do i=bfrag(3,j),bfrag(4,j)
-          mask_phi(i)=0
-          mask_theta(i)=0
-          iff(i)=1
-         enddo
-        else
-         do i=bfrag(4,j),bfrag(3,j)
-          mask_phi(i)=0
-          mask_theta(i)=0
-          iff(i)=1
-         enddo
-        endif
-       enddo
-       do j=1,nhfrag
-        do i=hfrag(1,j),hfrag(2,j)
-         mask_phi(i)=0
-         mask_theta(i)=0
-         iff(i)=1
-        enddo
-       enddo
-       mask_r=.true.
-
-
-
-       nhpb0=nhpb
-c
-c store dist. constrains
-c
-       do i=1,nres-3                                                             
-         do j=i+3,nres                                                           
-           if ( iff(i).eq.1.and.iff(j).eq.1 ) then
-            nhpb=nhpb+1                                                           
-            ihpb(nhpb)=i                                                          
-            jhpb(nhpb)=j                                                          
-            forcon(nhpb)=0.1                                                     
-            dhpb(nhpb)=DIST(i,j)
-           endif
-         enddo                                                                   
-       enddo                                    
-       call hpb_partition
-
-       if (debug) then
-        call chainbuild
-        call write_pdb(100+in_pdb,'input reg. structure',0d0)
-       endif
-       
-
-       ipot0=ipot
-       maxmin0=maxmin
-       maxfun0=maxfun
-       wstrain0=wstrain
-       wang0=wang
-c
-c      run soft pot. optimization 
-c
-       ipot=6
-       wang=3.0
-       maxmin=2000
-       maxfun=4000
-       call geom_to_var(nvar,var)
-#ifdef MPI
-       time0=MPI_WTIME()
-#else
-       time0=tcpu()
-#endif
-       call minimize(etot,var,iretcode,nfun)                               
-
-       write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun   
-#ifdef MPI
-       time1=MPI_WTIME()
-#else
-       time1=tcpu()
-#endif
-       write (iout,'(a,f6.2,f8.2,a)')'  Time for soft min.',time1-time0,
-     &         nfun/(time1-time0),' SOFT eval/s'
-       if (debug) then
-         call var_to_geom(nvar,var)
-         call chainbuild
-         call write_pdb(300+in_pdb,'soft structure',etot)
-       endif
-c
-c      run full UNRES optimization with constrains and frozen 2D
-c      the same variables as soft pot. optimizatio
-c
-       ipot=ipot0
-       wang=wang0
-       maxmin=maxmin0
-       maxfun=maxfun0
-#ifdef MPI
-       time0=MPI_WTIME()
-#else
-       time0=tcpu()
-#endif
-       call minimize(etot,var,iretcode,nfun)
-       write(iout,*)'SUMSL MASK DIST return code is',iretcode,
-     &                          ' eval ',nfun
-       ieval=nfun
-#ifdef MPI
-       time1=MPI_WTIME()
-#else
-       time1=tcpu()
-#endif
-       write (iout,'(a,f6.2,f8.2,a)') 
-     &        '  Time for mask dist min.',time1-time0,
-     &         nfun/(time1-time0),'  eval/s'
-       if (debug) then
-         call var_to_geom(nvar,var)
-         call chainbuild
-         call write_pdb(400+in_pdb,'mask & dist',etot)
-       endif
-c
-c      switch off constrains and 
-c      run full UNRES optimization with frozen 2D 
-c
-
-c
-c      reset constrains
-c
-       nhpb_c=nhpb
-       nhpb=nhpb0                                                                  
-       link_start=1                                                            
-       link_end=nhpb     
-       wstrain=wstrain0
-
-#ifdef MPI
-       time0=MPI_WTIME()
-#else
-       time0=tcpu()
-#endif
-       call minimize(etot,var,iretcode,nfun)
-       write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
-       ieval=ieval+nfun
-#ifdef MPI
-       time1=MPI_WTIME()
-#else
-       time1=tcpu()
-#endif
-       write (iout,'(a,f6.2,f8.2,a)')'  Time for mask min.',time1-time0,
-     &         nfun/(time1-time0),'  eval/s'
-
-
-       if (debug) then
-        call var_to_geom(nvar,var)
-        call chainbuild
-        call write_pdb(500+in_pdb,'mask 2d frozen',etot)
-       endif
-
-       mask_r=.false.
-
-
-c
-c      run full UNRES optimization with constrains and NO frozen 2D
-c
-
-       nhpb=nhpb_c                                                                  
-       link_start=1                                                            
-       link_end=nhpb     
-       maxfun=maxfun0/5
-
-       do ico=1,5
-
-       wstrain=wstrain0/ico
-#ifdef MPI
-       time0=MPI_WTIME()
-#else
-       time0=tcpu()
-#endif
-       call minimize(etot,var,iretcode,nfun)
-       write(iout,'(a10,f6.3,a14,i3,a6,i5)')
-     &   ' SUMSL DIST',wstrain,' return code is',iretcode,
-     &                          ' eval ',nfun
-       ieval=nfun
-#ifdef MPI
-       time1=MPI_WTIME()
-#else
-       time1=tcpu()
-#endif
-       write (iout,'(a,f6.2,f8.2,a)') 
-     &        '  Time for dist min.',time1-time0,
-     &         nfun/(time1-time0),'  eval/s'
-       if (debug) then
-         call var_to_geom(nvar,var)
-         call chainbuild
-         call write_pdb(600+in_pdb+ico,'dist cons',etot)
-       endif
-
-       enddo
-c
-       nhpb=nhpb0                                                                  
-       link_start=1                                                            
-       link_end=nhpb     
-       wstrain=wstrain0
-       maxfun=maxfun0
-
-
-c
-      if (minim) then
-#ifdef MPI
-       time0=MPI_WTIME()
-#else
-       time0=tcpu()
-#endif
-       call minimize(etot,var,iretcode,nfun)
-       write(iout,*)'------------------------------------------------'
-       write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
-     &  '+ DIST eval',ieval
-#ifdef MPI      
-       time1=MPI_WTIME()
-#else
-       time1=tcpu()
-#endif
-       write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
-     &         nfun/(time1-time0),' eval/s'
-
-
-       call var_to_geom(nvar,var)
-       call chainbuild        
-       call write_pdb(999,'full min',etot)
-      endif
-       
-      return
-      end
-
-