X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_Eshel%2FSRC-SURPLUS%2Ftest.F;fp=source%2Funres%2Fsrc_Eshel%2FSRC-SURPLUS%2Ftest.F;h=0140ee50675d9b8491e527ac5ef9aa6feb379eee;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_Eshel/SRC-SURPLUS/test.F b/source/unres/src_Eshel/SRC-SURPLUS/test.F new file mode 100644 index 0000000..0140ee5 --- /dev/null +++ b/source/unres/src_Eshel/SRC-SURPLUS/test.F @@ -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 + +