--- /dev/null
+ 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
+
+