Merge branch 'lipid' of mmka.chem.univ.gda.pl:unres into lipid
[unres.git] / source / unres / src_MD-NEWSC / test.F
diff --git a/source/unres/src_MD-NEWSC/test.F b/source/unres/src_MD-NEWSC/test.F
new file mode 100644 (file)
index 0000000..0140ee5
--- /dev/null
@@ -0,0 +1,863 @@
+      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
+
+