X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Ftest.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Ftest.F;h=4c7a728b010e4ca832ef15eecefc334dc4f2a969;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/test.F b/source/unres/src_MD-M-newcorr/test.F new file mode 100644 index 0000000..4c7a728 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/test.F @@ -0,0 +1,2707 @@ + subroutine test + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + 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_n16 + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + 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 jdata(5) + logical debug + debug=.true. + +c + call geom_to_var(nvar,var1) + call chainbuild + call etotal(energy(0)) + etot=energy(0) + write(iout,*) nnt,nct,etot + call write_pdb(1,'first structure',etot) + call secondary2(.true.) + + do i=1,4 + jdata(i)=bfrag(i,2) + enddo + + DO ij=1,4 + ieval=0 + jdata(5)=ij + call var_to_geom(nvar,var1) + write(iout,*) 'N16 test',(jdata(i),i=1,5) + call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5) + & ,ieval,ij) + call geom_to_var(nvar,var) + + if (minim) then + time0=MPI_WTIME() + call minimize(etot,var,iretcode,nfun) + write(iout,*)'------------------------------------------------' + write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, + & '+ DIST eval',ieval + + time1=MPI_WTIME() + 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(ij*100+99,'full min',etot) + endif + + + ENDDO + + 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 test11 + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.GEO' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.VAR' + include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' + include 'COMMON.FFIELD' + include 'COMMON.MINIM' +c + include 'COMMON.DISTFIT' + integer if(20,maxres),nif,ifa(20) + integer ibc(0:maxres,0:maxres),istrand(20) + integer ibd(maxres),ifb(10,2),nifb,lifb(10),lifb0 + integer itmp(20,maxres) + double precision time0,time1 + double precision energy(0:n_ene),ee + double precision varia(maxvar),vorg(maxvar) +c + logical debug,ltest,usedbfrag(maxres/3) + character*50 linia +c + integer betasheet(maxres),ibetasheet(maxres),nbetasheet + integer bstrand(maxres/3,6),nbstrand + +c------------------------ + + debug=.true. +c------------------------ + nbstrand=0 + nbetasheet=0 + do i=1,nres + betasheet(i)=0 + ibetasheet(i)=0 + enddo + call geom_to_var(nvar,vorg) + call secondary2(debug) + + if (nbfrag.le.1) return + + do i=1,nbfrag + usedbfrag(i)=.false. + enddo + + + nbetasheet=nbetasheet+1 + nbstrand=2 + bstrand(1,1)=bfrag(1,1) + bstrand(1,2)=bfrag(2,1) + bstrand(1,3)=nbetasheet + bstrand(1,4)=1 + bstrand(1,5)=bfrag(1,1) + bstrand(1,6)=bfrag(2,1) + do i=bfrag(1,1),bfrag(2,1) + betasheet(i)=nbetasheet + ibetasheet(i)=1 + enddo +c + bstrand(2,1)=bfrag(3,1) + bstrand(2,2)=bfrag(4,1) + bstrand(2,3)=nbetasheet + bstrand(2,5)=bfrag(3,1) + bstrand(2,6)=bfrag(4,1) + + if (bfrag(3,1).le.bfrag(4,1)) then + bstrand(2,4)=2 + do i=bfrag(3,1),bfrag(4,1) + betasheet(i)=nbetasheet + ibetasheet(i)=2 + enddo + else + bstrand(2,4)=-2 + do i=bfrag(4,1),bfrag(3,1) + betasheet(i)=nbetasheet + ibetasheet(i)=2 + enddo + endif + + iused_nbfrag=1 + + do while (iused_nbfrag.ne.nbfrag) + + do j=2,nbfrag + + IF (.not.usedbfrag(j)) THEN + + write (*,*) j,(bfrag(i,j),i=1,4) + do jk=6,1,-1 + write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand) + enddo + write (*,*) '------------------' + + + if (bfrag(3,j).le.bfrag(4,j)) then + do i=bfrag(3,j),bfrag(4,j) + if(betasheet(i).eq.nbetasheet) then + in=ibetasheet(i) + do k=bfrag(3,j),bfrag(4,j) + betasheet(k)=nbetasheet + ibetasheet(k)=in + enddo + nbstrand=nbstrand+1 + usedbfrag(j)=.true. + iused_nbfrag=iused_nbfrag+1 + do k=bfrag(1,j),bfrag(2,j) + betasheet(k)=nbetasheet + ibetasheet(k)=nbstrand + enddo + if (bstrand(in,4).lt.0) then + bstrand(nbstrand,1)=bfrag(2,j) + bstrand(nbstrand,2)=bfrag(1,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,4)=-nbstrand + bstrand(nbstrand,5)=bstrand(nbstrand,1) + bstrand(nbstrand,6)=bstrand(nbstrand,2) + if(bstrand(in,1).lt.bfrag(4,j)) then + call bgrow(bstrand,nbstrand,in,1,bfrag(4,j)) + else + bstrand(nbstrand,5)=bstrand(nbstrand,5)+ + & (bstrand(in,5)-bfrag(4,j)) + endif + if(bstrand(in,2).gt.bfrag(3,j)) then + call bgrow(bstrand,nbstrand,in,2,bfrag(3,j)) + else + bstrand(nbstrand,6)=bstrand(nbstrand,6)- + & (-bstrand(in,6)+bfrag(3,j)) + endif + else + bstrand(nbstrand,1)=bfrag(1,j) + bstrand(nbstrand,2)=bfrag(2,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,4)=nbstrand + bstrand(nbstrand,5)=bstrand(nbstrand,1) + bstrand(nbstrand,6)=bstrand(nbstrand,2) + if(bstrand(in,1).gt.bfrag(3,j)) then + call bgrow(bstrand,nbstrand,in,1,bfrag(3,j)) + else + bstrand(nbstrand,5)=bstrand(nbstrand,5)- + & (-bstrand(in,5)+bfrag(3,j)) + endif + if(bstrand(in,2).lt.bfrag(4,j)) then + call bgrow(bstrand,nbstrand,in,2,bfrag(4,j)) + else + bstrand(nbstrand,6)=bstrand(nbstrand,6)+ + & (bstrand(in,6)-bfrag(4,j)) + endif + endif + goto 11 + endif + if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then + in=ibetasheet(bfrag(1,j)+i-bfrag(3,j)) + do k=bfrag(1,j),bfrag(2,j) + betasheet(k)=nbetasheet + ibetasheet(k)=in + enddo + nbstrand=nbstrand+1 + usedbfrag(j)=.true. + iused_nbfrag=iused_nbfrag+1 + do k=bfrag(3,1),bfrag(4,1) + betasheet(k)=nbetasheet + ibetasheet(k)=nbstrand + enddo + if (bstrand(in,4).lt.0) then + bstrand(nbstrand,1)=bfrag(4,j) + bstrand(nbstrand,2)=bfrag(3,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,4)=-nbstrand + bstrand(nbstrand,5)=bstrand(nbstrand,1) + bstrand(nbstrand,6)=bstrand(nbstrand,2) + if(bstrand(in,1).lt.bfrag(2,j)) then + call bgrow(bstrand,nbstrand,in,1,bfrag(2,j)) + else + bstrand(nbstrand,5)=bstrand(nbstrand,5)+ + & (bstrand(in,5)-bfrag(2,j)) + endif + if(bstrand(in,2).gt.bfrag(1,j)) then + call bgrow(bstrand,nbstrand,in,2,bfrag(1,j)) + else + bstrand(nbstrand,6)=bstrand(nbstrand,6)- + & (-bstrand(in,6)+bfrag(1,j)) + endif + else + bstrand(nbstrand,1)=bfrag(3,j) + bstrand(nbstrand,2)=bfrag(4,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,4)=nbstrand + bstrand(nbstrand,5)=bstrand(nbstrand,1) + bstrand(nbstrand,6)=bstrand(nbstrand,2) + if(bstrand(in,1).gt.bfrag(1,j)) then + call bgrow(bstrand,nbstrand,in,1,bfrag(1,j)) + else + bstrand(nbstrand,5)=bstrand(nbstrand,5)- + & (-bstrand(in,5)+bfrag(1,j)) + endif + if(bstrand(in,2).lt.bfrag(2,j)) then + call bgrow(bstrand,nbstrand,in,2,bfrag(2,j)) + else + bstrand(nbstrand,6)=bstrand(nbstrand,6)+ + & (bstrand(in,6)-bfrag(2,j)) + endif + endif + goto 11 + endif + enddo + else + do i=bfrag(4,j),bfrag(3,j) + if(betasheet(i).eq.nbetasheet) then + in=ibetasheet(i) + do k=bfrag(4,j),bfrag(3,j) + betasheet(k)=nbetasheet + ibetasheet(k)=in + enddo + nbstrand=nbstrand+1 + usedbfrag(j)=.true. + iused_nbfrag=iused_nbfrag+1 + do k=bfrag(1,j),bfrag(2,j) + betasheet(k)=nbetasheet + ibetasheet(k)=nbstrand + enddo + if (bstrand(in,4).lt.0) then + bstrand(nbstrand,1)=bfrag(1,j) + bstrand(nbstrand,2)=bfrag(2,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,4)=nbstrand + bstrand(nbstrand,5)=bstrand(nbstrand,1) + bstrand(nbstrand,6)=bstrand(nbstrand,2) + if(bstrand(in,1).lt.bfrag(3,j)) then + call bgrow(bstrand,nbstrand,in,1,bfrag(3,j)) + else + bstrand(nbstrand,5)=bstrand(nbstrand,5)- + & (bstrand(in,5)-bfrag(3,j)) + endif + if(bstrand(in,2).gt.bfrag(4,j)) then + call bgrow(bstrand,nbstrand,in,2,bfrag(4,j)) + else + bstrand(nbstrand,6)=bstrand(nbstrand,6)+ + & (-bstrand(in,6)+bfrag(4,j)) + endif + else + bstrand(nbstrand,1)=bfrag(2,j) + bstrand(nbstrand,2)=bfrag(1,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,4)=-nbstrand + bstrand(nbstrand,5)=bstrand(nbstrand,1) + bstrand(nbstrand,6)=bstrand(nbstrand,2) + if(bstrand(in,1).gt.bfrag(4,j)) then + call bgrow(bstrand,nbstrand,in,1,bfrag(4,j)) + else + bstrand(nbstrand,5)=bstrand(nbstrand,5)+ + & (-bstrand(in,5)+bfrag(4,j)) + endif + if(bstrand(in,2).lt.bfrag(3,j)) then + call bgrow(bstrand,nbstrand,in,2,bfrag(3,j)) + else + bstrand(nbstrand,6)=bstrand(nbstrand,6)- + & (bstrand(in,6)-bfrag(3,j)) + endif + endif + goto 11 + endif + if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then + in=ibetasheet(bfrag(2,j)-i+bfrag(4,j)) + do k=bfrag(1,j),bfrag(2,j) + betasheet(k)=nbetasheet + ibetasheet(k)=in + enddo + nbstrand=nbstrand+1 + usedbfrag(j)=.true. + iused_nbfrag=iused_nbfrag+1 + do k=bfrag(4,j),bfrag(3,j) + betasheet(k)=nbetasheet + ibetasheet(k)=nbstrand + enddo + if (bstrand(in,4).lt.0) then + bstrand(nbstrand,1)=bfrag(4,j) + bstrand(nbstrand,2)=bfrag(3,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,4)=nbstrand + bstrand(nbstrand,5)=bstrand(nbstrand,1) + bstrand(nbstrand,6)=bstrand(nbstrand,2) + if(bstrand(in,1).lt.bfrag(2,j)) then + call bgrow(bstrand,nbstrand,in,1,bfrag(2,j)) + else + bstrand(nbstrand,5)=bstrand(nbstrand,5)- + & (bstrand(in,5)-bfrag(2,j)) + endif + if(bstrand(in,2).gt.bfrag(1,j)) then + call bgrow(bstrand,nbstrand,in,2,bfrag(1,j)) + else + bstrand(nbstrand,6)=bstrand(nbstrand,6)+ + & (-bstrand(in,6)+bfrag(1,j)) + endif + else + bstrand(nbstrand,1)=bfrag(3,j) + bstrand(nbstrand,2)=bfrag(4,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,4)=-nbstrand + bstrand(nbstrand,5)=bstrand(nbstrand,1) + bstrand(nbstrand,6)=bstrand(nbstrand,2) + if(bstrand(in,1).gt.bfrag(1,j)) then + call bgrow(bstrand,nbstrand,in,1,bfrag(1,j)) + else + bstrand(nbstrand,5)=bstrand(nbstrand,5)+ + & (-bstrand(in,5)+bfrag(1,j)) + endif + if(bstrand(in,2).lt.bfrag(2,j)) then + call bgrow(bstrand,nbstrand,in,2,bfrag(2,j)) + else + bstrand(nbstrand,6)=bstrand(nbstrand,6)- + & (bstrand(in,6)-bfrag(2,j)) + endif + endif + goto 11 + endif + enddo + endif + + + + ENDIF + enddo + + j=2 + do while (usedbfrag(j)) + j=j+1 + enddo + + nbstrand=nbstrand+1 + nbetasheet=nbetasheet+1 + bstrand(nbstrand,1)=bfrag(1,j) + bstrand(nbstrand,2)=bfrag(2,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,5)=bfrag(1,j) + bstrand(nbstrand,6)=bfrag(2,j) + + bstrand(nbstrand,4)=nbstrand + do i=bfrag(1,j),bfrag(2,j) + betasheet(i)=nbetasheet + ibetasheet(i)=nbstrand + enddo +c + nbstrand=nbstrand+1 + bstrand(nbstrand,1)=bfrag(3,j) + bstrand(nbstrand,2)=bfrag(4,j) + bstrand(nbstrand,3)=nbetasheet + bstrand(nbstrand,5)=bfrag(3,j) + bstrand(nbstrand,6)=bfrag(4,j) + + if (bfrag(3,j).le.bfrag(4,j)) then + bstrand(nbstrand,4)=nbstrand + do i=bfrag(3,j),bfrag(4,j) + betasheet(i)=nbetasheet + ibetasheet(i)=nbstrand + enddo + else + bstrand(nbstrand,4)=-nbstrand + do i=bfrag(4,j),bfrag(3,j) + betasheet(i)=nbetasheet + ibetasheet(i)=nbstrand + enddo + endif + + iused_nbfrag=iused_nbfrag+1 + usedbfrag(j)=.true. + + + 11 continue + do jk=6,1,-1 + write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand) + enddo + + + enddo + + do i=1,nres + if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i) + enddo + write(*,*) + do j=6,1,-1 + write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand) + enddo + +c------------------------ + nifb=0 + do i=1,nbstrand + do j=i+1,nbstrand + if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or. + & iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then + nifb=nifb+1 + ifb(nifb,1)=bstrand(i,4) + ifb(nifb,2)=bstrand(j,4) + endif + enddo + enddo + + write(*,*) + do i=1,nifb + write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2) + enddo + + do i=1,nbstrand + ifa(i)=bstrand(i,4) + enddo + write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand) + + nif=iabs(bstrand(1,6)-bstrand(1,5))+1 + do j=2,nbstrand + if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif) + & nif=iabs(bstrand(j,6)-bstrand(j,5))+1 + enddo + + write(*,*) nif + do i=1,nif + do j=1,nbstrand + if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6)) + if (if(j,i).gt.0) then + if(betasheet(if(j,i)).eq.0 .or. + & ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0 + else + if(j,i)=0 + endif + enddo + write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand) + enddo + +c read (inp,*) (ifa(i),i=1,4) +c do i=1,nres +c read (inp,*,err=20,end=20) (if(j,i),j=1,4) +c enddo +c 20 nif=i-1 + stop +c------------------------ + + isa=4 + is=2*isa-1 + iconf=0 +cccccccccccccccccccccccccccccccccc + DO ig=1,is**isa-1 +cccccccccccccccccccccccccccccccccc + + ii=ig + do j=1,is + istrand(is-j+1)=int(ii/is**(is-j)) + ii=ii-istrand(is-j+1)*is**(is-j) + enddo + ltest=.true. + do k=1,isa + istrand(k)=istrand(k)+1 + if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1 + enddo + do k=1,isa + do l=1,isa + if(istrand(k).eq.istrand(l).and.k.ne.l.or. + & istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false. + enddo + enddo + + lifb0=1 + do m=1,nifb + lifb(m)=0 + do k=1,isa-1 + if( + & ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. + & ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. + & -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. + & -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) + & lifb(m)=1 + enddo + lifb0=lifb0*lifb(m) + enddo + + if (mod(isa,2).eq.0) then + do k=isa/2+1,isa + if (istrand(k).eq.1) ltest=.false. + enddo + else + do k=(isa+1)/2+1,isa + if (istrand(k).eq.1) ltest=.false. + enddo + endif + + IF (ltest.and.lifb0.eq.1) THEN + iconf=iconf+1 + + call var_to_geom(nvar,vorg) + + write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa) + write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa) + write (linia,'(10i3)') (istrand(k),k=1,isa) + + do i=1,nres + do j=1,nres + ibc(i,j)=0 + enddo + enddo + + + do i=1,4 + if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then + do j=1,nif + itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j) + enddo + else + do j=1,nif + itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1) + enddo + endif + enddo + + do i=1,nif + write(*,*) (itmp(j,i),j=1,4) + enddo + + do i=1,nif +c ifa(1),ifa(2),ifa(3),ifa(4) +c if(1,i),if(2,i),if(3,i),if(4,i) + do k=1,isa-1 + ltest=.false. + do m=1,nifb + if( + & ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or. + & ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or. + & -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or. + & -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1)) + & then + ltest=.true. + goto 110 + endif + enddo + 110 continue + if (ltest) then + ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1 + else + ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2 + endif +c + if (k.lt.3) + & ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3 + if (k.lt.2) + & ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4 + enddo + enddo +c------------------------ + +c +c freeze sec.elements +c + do i=1,nres + mask(i)=1 + mask_phi(i)=1 + mask_theta(i)=1 + mask_side(i)=1 + enddo + + do j=1,nbfrag + do i=bfrag(1,j),bfrag(2,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + if (bfrag(3,j).le.bfrag(4,j)) then + do i=bfrag(3,j),bfrag(4,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + else + do i=bfrag(4,j),bfrag(3,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + endif + enddo + do j=1,nhfrag + do i=hfrag(1,j),hfrag(2,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + enddo + mask_r=.true. + +c------------------------ +c generate constrains +c + nhpb0=nhpb + call chainbuild + ind=0 + do i=1,nres-3 + do j=i+3,nres + ind=ind+1 + if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then + d0(ind)=DIST(i,j) + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then + d0(ind)=5.0 + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then + d0(ind)=11.0 + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then + d0(ind)=16.0 + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else if ( ibc(i,j).gt.0 ) then + d0(ind)=DIST(i,ibc(i,j)) + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else if ( ibc(j,i).gt.0 ) then + d0(ind)=DIST(ibc(j,i),j) + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else + w(ind)=0.0 + endif + dd(ind)=d0(ind) + enddo + enddo + call hpb_partition +cd-------------------------- + + write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)), + & ibc(jhpb(i),ihpb(i)),' --', + & ihpb(i),jhpb(i),dhpb(i),i=1,nhpb) + +cd nhpb=0 +cd goto 901 +c +c + call contact_cp_min(varia,ifun,iconf,linia,debug) + if (minim) then + time0=MPI_WTIME() + call minimize(etot,varia,iretcode,nfun) + write(iout,*)'------------------------------------------------' + write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, + & '+ DIST eval',ifun + + time1=MPI_WTIME() + write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0, + & nfun/(time1-time0),' eval/s' + + write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa) + call var_to_geom(nvar,varia) + call chainbuild + call write_pdb(900+iconf,linia,etot) + endif + + call etotal(energy(0)) + etot=energy(0) + call enerprint(energy(0)) +cd call intout +cd call briefout(0,etot) +cd call secondary2(.true.) + + 901 CONTINUE +ctest return +cccccccccccccccccccccccccccccccccccc + ENDIF + ENDDO +cccccccccccccccccccccccccccccccccccc + + return + 10 write (iout,'(a)') 'Error reading test structure.' + return + end +c-------------------------------------------------------- + + subroutine test3 + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.GEO' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.VAR' + include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' + include 'COMMON.FFIELD' + include 'COMMON.MINIM' +c + include 'COMMON.DISTFIT' + integer if(3,maxres),nif + integer ibc(maxres,maxres),istrand(20) + integer ibd(maxres),ifb(10,2),nifb,lifb(10),lifb0 + double precision time0,time1 + double precision energy(0:n_ene),ee + double precision varia(maxvar) +c + logical debug,ltest + character*50 linia +c + do i=1,nres + read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i) + enddo + 20 nif=i-1 + write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i), + & i=1,nif) + + +c------------------------ + call secondary2(debug) +c------------------------ + do i=1,nres + do j=1,nres + ibc(i,j)=0 + enddo + enddo + +c +c freeze sec.elements and store indexes for beta constrains +c + do i=1,nres + mask(i)=1 + mask_phi(i)=1 + mask_theta(i)=1 + mask_side(i)=1 + enddo + + do j=1,nbfrag + do i=bfrag(1,j),bfrag(2,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + if (bfrag(3,j).le.bfrag(4,j)) then + do i=bfrag(3,j),bfrag(4,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1 + enddo + else + do i=bfrag(4,j),bfrag(3,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1 + enddo + endif + enddo + do j=1,nhfrag + do i=hfrag(1,j),hfrag(2,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + enddo + mask_r=.true. + + +c ---------------- test -------------- + do i=1,nif + if (ibc(if(1,i),if(2,i)).eq.-1) then + ibc(if(1,i),if(2,i))=if(3,i) + ibc(if(1,i),if(3,i))=if(2,i) + else if (ibc(if(2,i),if(1,i)).eq.-1) then + ibc(if(2,i),if(1,i))=0 + ibc(if(1,i),if(2,i))=if(3,i) + ibc(if(1,i),if(3,i))=if(2,i) + else + ibc(if(1,i),if(2,i))=if(3,i) + ibc(if(1,i),if(3,i))=if(2,i) + endif + enddo + + do i=1,nres + do j=1,nres + if (ibc(i,j).ne.0) write(*,'(3i5)') i,j,ibc(i,j) + enddo + enddo +c------------------------ + call chainbuild + ind=0 + do i=1,nres-3 + do j=i+3,nres + ind=ind+1 + if ( ibc(i,j).eq.-1 ) then + d0(ind)=DIST(i,j) + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else if ( ibc(i,j).gt.0 ) then + d0(ind)=DIST(i,ibc(i,j)) + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else if ( ibc(j,i).gt.0 ) then + d0(ind)=DIST(ibc(j,i),j) + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else + w(ind)=0.0 + endif + enddo + enddo + call hpb_partition + +cd-------------------------- + write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)), + & ibc(jhpb(i),ihpb(i)),' --', + & ihpb(i),jhpb(i),dhpb(i),i=1,nhpb) + + + linia='dist' + debug=.true. + in_pdb=7 +c + call contact_cp_min(varia,ieval,in_pdb,linia,debug) + if (minim) then + time0=MPI_WTIME() + call minimize(etot,varia,iretcode,nfun) + write(iout,*)'------------------------------------------------' + write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, + & '+ DIST eval',ieval + + time1=MPI_WTIME() + write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0, + & nfun/(time1-time0),' eval/s' + + + call var_to_geom(nvar,varia) + call chainbuild + call write_pdb(999,'full min',etot) + endif + + call etotal(energy(0)) + etot=energy(0) + call enerprint(energy(0)) + call intout + call briefout(0,etot) + call secondary2(.true.) + + return + 10 write (iout,'(a)') 'Error reading test structure.' + return + end + + + + + subroutine test__ + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.GEO' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.VAR' + include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' + include 'COMMON.FFIELD' + include 'COMMON.MINIM' +c + include 'COMMON.DISTFIT' + integer if(2,2),ind + integer iff(maxres) + double precision time0,time1 + double precision energy(0:n_ene),ee + double precision theta2(maxres),phi2(maxres),alph2(maxres), + & omeg2(maxres), + & theta1(maxres),phi1(maxres),alph1(maxres), + & omeg1(maxres) + double precision varia(maxvar),varia2(maxvar) +c + + + read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2) + write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2) + read (inp,*,err=10,end=10) (theta2(i),i=3,nres) + read (inp,*,err=10,end=10) (phi2(i),i=4,nres) + read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1) + read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1) + do i=1,nres + theta2(i)=deg2rad*theta2(i) + phi2(i)=deg2rad*phi2(i) + alph2(i)=deg2rad*alph2(i) + omeg2(i)=deg2rad*omeg2(i) + enddo + do i=1,nres + theta1(i)=theta(i) + phi1(i)=phi(i) + alph1(i)=alph(i) + omeg1(i)=omeg(i) + enddo + + do i=1,nres + mask(i)=1 + enddo + + +c------------------------ + do i=1,nres + iff(i)=0 + enddo + do j=1,2 + do i=if(j,1),if(j,2) + iff(i)=1 + enddo + enddo + + call chainbuild + call geom_to_var(nvar,varia) + call write_pdb(1,'first structure',0d0) + + call secondary(.true.) + + call secondary2(.true.) + + do j=1,nbfrag + if ( (bfrag(3,j).lt.bfrag(4,j) .or. + & bfrag(4,j)-bfrag(2,j).gt.4) .and. + & bfrag(2,j)-bfrag(1,j).gt.3 ) then + nn=nn+1 + + if (bfrag(3,j).lt.bfrag(4,j)) then + write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') + & "select",bfrag(1,j)-1,"-",bfrag(2,j)-1 + & ,",",bfrag(3,j)-1,"-",bfrag(4,j)-1 + else + write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') + & "select",bfrag(1,j)-1,"-",bfrag(2,j)-1 + & ,",",bfrag(4,j)-1,"-",bfrag(3,j)-1 + endif + endif + enddo + + do i=1,nres + theta(i)=theta2(i) + phi(i)=phi2(i) + alph(i)=alph2(i) + omeg(i)=omeg2(i) + enddo + + call chainbuild + call geom_to_var(nvar,varia2) + call write_pdb(2,'second structure',0d0) + + + +c------------------------------------------------------- + + ifun=-1 + call contact_cp(varia,varia2,iff,ifun,7) + if (minim) then + time0=MPI_WTIME() + call minimize(etot,varia,iretcode,nfun) + write(iout,*)'------------------------------------------------' + write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, + & '+ DIST eval',ifun + + time1=MPI_WTIME() + write (iout,'(a,f6.2,f8.2,a)')' Time for full min.',time1-time0, + & nfun/(time1-time0),' eval/s' + + + call var_to_geom(nvar,varia) + call chainbuild + call write_pdb(999,'full min',etot) + endif + + call etotal(energy(0)) + etot=energy(0) + call enerprint(energy(0)) + call intout + call briefout(0,etot) + + return + 10 write (iout,'(a)') 'Error reading test structure.' + return + end + +c------------------------------------------------- +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 contact_cp2(var,var2,iff,ieval,in_pdb) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.SBRIDGE' + include 'COMMON.FFIELD' + include 'COMMON.IOUNITS' + include 'COMMON.DISTFIT' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.MINIM' + + character*50 linia + integer nf,ij(4) + double precision var(maxvar),var2(maxvar) + double precision time0,time1 + integer iff(maxres),ieval + double precision theta1(maxres),phi1(maxres),alph1(maxres), + & omeg1(maxres) + + + call var_to_geom(nvar,var) + call chainbuild + nhpb0=nhpb + ind=0 + do i=1,nres-3 + do j=i+3,nres + ind=ind+1 + if ( iff(i).eq.1.and.iff(j).eq.1 ) then + d0(ind)=DIST(i,j) + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else + w(ind)=0.0 + endif + enddo + enddo + call hpb_partition + + do i=1,nres + theta1(i)=theta(i) + phi1(i)=phi(i) + alph1(i)=alph(i) + omeg1(i)=omeg(i) + enddo + + call var_to_geom(nvar,var2) + + do i=1,nres + if ( iff(i).eq.1 ) then + theta(i)=theta1(i) + phi(i)=phi1(i) + alph(i)=alph1(i) + omeg(i)=omeg1(i) + endif + enddo + + call chainbuild +cd call write_pdb(3,'combined structure',0d0) +cd time0=MPI_WTIME() + + NX=NRES-3 + NY=((NRES-4)*(NRES-5))/2 + call distfit(.true.,200) + +cd time1=MPI_WTIME() +cd write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec' + + ipot0=ipot + maxmin0=maxmin + maxfun0=maxfun + wstrain0=wstrain + + ipot=6 + maxmin=2000 + maxfun=5000 + call geom_to_var(nvar,var) +cd time0=MPI_WTIME() + call minimize(etot,var,iretcode,nfun) + write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun + +cd time1=MPI_WTIME() +cd write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0, +cd & nfun/(time1-time0),' SOFT eval/s' + call var_to_geom(nvar,var) + call chainbuild + + + iwsk=0 + nf=0 + if (iff(1).eq.1) then + iwsk=1 + nf=nf+1 + ij(nf)=0 + endif + do i=2,nres + if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then + iwsk=1 + nf=nf+1 + ij(nf)=i + endif + if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then + iwsk=0 + nf=nf+1 + ij(nf)=i-1 + endif + enddo + if (iff(nres).eq.1) then + nf=nf+1 + ij(nf)=nres + endif + + +cd write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') +cd & "select",ij(1),"-",ij(2), +cd & ",",ij(3),"-",ij(4) +cd call write_pdb(in_pdb,linia,etot) + + + ipot=ipot0 + maxmin=maxmin0 + maxfun=maxfun0 +cd time0=MPI_WTIME() + call minimize(etot,var,iretcode,nfun) +cd write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun + ieval=nfun + +cd time1=MPI_WTIME() +cd write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0, +cd & nfun/(time1-time0),' eval/s' +cd call var_to_geom(nvar,var) +cd call chainbuild +cd call write_pdb(6,'dist structure',etot) + + + nhpb= nhpb0 + link_start=1 + link_end=nhpb + wstrain=wstrain0 + + return + end +c----------------------------------------------------------- + subroutine contact_cp(var,var2,iff,ieval,in_pdb) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.SBRIDGE' + include 'COMMON.FFIELD' + include 'COMMON.IOUNITS' + include 'COMMON.DISTFIT' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.MINIM' + + character*50 linia + integer nf,ij(4) + double precision energy(0:n_ene) + double precision var(maxvar),var2(maxvar) + double precision time0,time1 + integer iff(maxres),ieval + double precision theta1(maxres),phi1(maxres),alph1(maxres), + & omeg1(maxres) + logical debug + + debug=.false. +c debug=.true. + if (ieval.eq.-1) debug=.true. + + +c +c store selected dist. constrains from 1st structure +c +#ifdef OSF +c Intercept NaNs in the coordinates +c write(iout,*) (var(i),i=1,nvar) + x_sum=0.D0 + do i=1,nvar + x_sum=x_sum+var(i) + enddo + if (x_sum.ne.x_sum) then + write(iout,*)" *** contact_cp : Found NaN in coordinates" + call flush(iout) + print *," *** contact_cp : Found NaN in coordinates" + return + endif +#endif + + + call var_to_geom(nvar,var) + call chainbuild + nhpb0=nhpb + ind=0 + do i=1,nres-3 + do j=i+3,nres + ind=ind+1 + if ( iff(i).eq.1.and.iff(j).eq.1 ) then + d0(ind)=DIST(i,j) + w(ind)=10.0 + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=10.0 + dhpb(nhpb)=d0(ind) + else + w(ind)=0.0 + endif + enddo + enddo + call hpb_partition + + do i=1,nres + theta1(i)=theta(i) + phi1(i)=phi(i) + alph1(i)=alph(i) + omeg1(i)=omeg(i) + enddo + +c +c freeze sec.elements from 2nd structure +c + do i=1,nres + mask_phi(i)=1 + mask_theta(i)=1 + mask_side(i)=1 + enddo + + call var_to_geom(nvar,var2) + call secondary2(debug) + do j=1,nbfrag + do i=bfrag(1,j),bfrag(2,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + if (bfrag(3,j).le.bfrag(4,j)) then + do i=bfrag(3,j),bfrag(4,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + else + do i=bfrag(4,j),bfrag(3,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + endif + enddo + do j=1,nhfrag + do i=hfrag(1,j),hfrag(2,j) + mask(i)=0 + mask_phi(i)=0 + mask_theta(i)=0 + enddo + enddo + mask_r=.true. + +c +c copy selected res from 1st to 2nd structure +c + + do i=1,nres + if ( iff(i).eq.1 ) then + theta(i)=theta1(i) + phi(i)=phi1(i) + alph(i)=alph1(i) + omeg(i)=omeg1(i) + endif + enddo + + if(debug) then +c +c prepare description in linia variable +c + iwsk=0 + nf=0 + if (iff(1).eq.1) then + iwsk=1 + nf=nf+1 + ij(nf)=1 + endif + do i=2,nres + if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then + iwsk=1 + nf=nf+1 + ij(nf)=i + endif + if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then + iwsk=0 + nf=nf+1 + ij(nf)=i-1 + endif + enddo + if (iff(nres).eq.1) then + nf=nf+1 + ij(nf)=nres + endif + + write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') + & "SELECT",ij(1)-1,"-",ij(2)-1, + & ",",ij(3)-1,"-",ij(4)-1 + + endif +c +c run optimization +c + call contact_cp_min(var,ieval,in_pdb,linia,debug) + + return + end + + subroutine contact_cp_min(var,ieval,in_pdb,linia,debug) +c +c input : theta,phi,alph,omeg,in_pdb,linia,debug +c output : var,ieval +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.SBRIDGE' + include 'COMMON.FFIELD' + include 'COMMON.IOUNITS' + include 'COMMON.DISTFIT' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.MINIM' + + character*50 linia + integer nf,ij(4) + double precision energy(0:n_ene) + double precision var(maxvar) + double precision time0,time1 + integer ieval,info(3) + logical debug,fail,check_var,reduce,change + + write(iout,'(a20,i6,a20)') + & '------------------',in_pdb,'-------------------' + + if (debug) then + call chainbuild + call write_pdb(1000+in_pdb,'combined structure',0d0) + time0=MPI_WTIME() + endif + +c +c run optimization of distances +c +c uses d0(),w() and mask() for frozen 2D +c +ctest--------------------------------------------- +ctest NX=NRES-3 +ctest NY=((NRES-4)*(NRES-5))/2 +ctest call distfit(debug,5000) + + do i=1,nres + mask_side(i)=0 + enddo + + ipot01=ipot + maxmin01=maxmin + maxfun01=maxfun +c wstrain01=wstrain + wsc01=wsc + wscp01=wscp + welec01=welec + wvdwpp01=wvdwpp +c wang01=wang + wscloc01=wscloc + wtor01=wtor + wtor_d01=wtor_d + + ipot=6 + maxmin=2000 + maxfun=4000 +c wstrain=1.0 + wsc=0.0 + wscp=0.0 + welec=0.0 + wvdwpp=0.0 +c wang=0.0 + wscloc=0.0 + wtor=0.0 + wtor_d=0.0 + + call geom_to_var(nvar,var) +cde change=reduce(var) + if (check_var(var,info)) then + write(iout,*) 'cp_min error in input' + print *,'cp_min error in input' + return + endif + +cd call etotal(energy(0)) +cd call enerprint(energy(0)) +cd call check_eint + + time0=MPI_WTIME() +cdtest call minimize(etot,var,iretcode,nfun) +cdtest write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun + time1=MPI_WTIME() + +cd call etotal(energy(0)) +cd call enerprint(energy(0)) +cd call check_eint + + do i=1,nres + mask_side(i)=1 + enddo + + ipot=ipot01 + maxmin=maxmin01 + maxfun=maxfun01 +c wstrain=wstrain01 + wsc=wsc01 + wscp=wscp01 + welec=welec01 + wvdwpp=wvdwpp01 +c wang=wang01 + wscloc=wscloc01 + wtor=wtor01 + wtor_d=wtor_d01 +ctest-------------------------------------------------- + + if(debug) then + time1=MPI_WTIME() + write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec' + call write_pdb(2000+in_pdb,'distfit structure',0d0) + endif + + + ipot0=ipot + maxmin0=maxmin + maxfun0=maxfun + wstrain0=wstrain +c +c run soft pot. optimization +c with constrains: +c nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition +c and frozen 2D: +c mask_phi(),mask_theta(),mask_side(),mask_r +c + ipot=6 + maxmin=2000 + maxfun=4000 + +cde change=reduce(var) +cde if (check_var(var,info)) write(iout,*) 'error before soft' + time0=MPI_WTIME() + call minimize(etot,var,iretcode,nfun) + + write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun + time1=MPI_WTIME() + 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(3000+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 + maxmin=maxmin0 + maxfun=maxfun0 +c +c check overlaps before calling full UNRES minim +c + call var_to_geom(nvar,var) + call chainbuild + call etotal(energy(0)) +#ifdef OSF + write(iout,*) 'N7 ',energy(0) + if (energy(0).ne.energy(0)) then + write(iout,*) 'N7 error - gives NaN',energy(0) + endif +#endif + ieval=1 + if (energy(1).eq.1.0d20) then + write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1) + call overlap_sc(fail) + if(.not.fail) then + call etotal(energy(0)) + ieval=ieval+1 + write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1) + else + mask_r=.false. + nhpb= nhpb0 + link_start=1 + link_end=nhpb + wstrain=wstrain0 + return + endif + endif + call flush(iout) +c +cdte time0=MPI_WTIME() +cde change=reduce(var) +cde if (check_var(var,info)) then +cde write(iout,*) 'error before mask dist' +cde call var_to_geom(nvar,var) +cde call chainbuild +cde call write_pdb(10000+in_pdb,'before mask dist',etot) +cde endif +cdte call minimize(etot,var,iretcode,nfun) +cdte write(iout,*)'SUMSL MASK DIST return code is',iretcode, +cdte & ' eval ',nfun +cdte ieval=ieval+nfun +cdte +cdte time1=MPI_WTIME() +cdte write (iout,'(a,f6.2,f8.2,a)') +cdte & ' Time for mask dist min.',time1-time0, +cdte & nfun/(time1-time0),' eval/s' +cdte call flush(iout) + if (debug) then + call var_to_geom(nvar,var) + call chainbuild + call write_pdb(4000+in_pdb,'mask dist',etot) + endif +c +c switch off freezing of 2D and +c run full UNRES optimization with constrains +c + mask_r=.false. + time0=MPI_WTIME() +cde change=reduce(var) +cde if (check_var(var,info)) then +cde write(iout,*) 'error before dist' +cde call var_to_geom(nvar,var) +cde call chainbuild +cde call write_pdb(11000+in_pdb,'before dist',etot) +cde endif + + call minimize(etot,var,iretcode,nfun) + +cde change=reduce(var) +cde if (check_var(var,info)) then +cde write(iout,*) 'error after dist',ico +cde call var_to_geom(nvar,var) +cde call chainbuild +cde call write_pdb(12000+in_pdb+ico*1000,'after dist',etot) +cde endif + write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun + ieval=ieval+nfun + + time1=MPI_WTIME() + write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0, + & nfun/(time1-time0),' eval/s' +cde call etotal(energy(0)) +cde write(iout,*) 'N7 after dist',energy(0) + call flush(iout) + + if (debug) then + call var_to_geom(nvar,var) + call chainbuild + call write_pdb(in_pdb,linia,etot) + endif +c +c reset constrains +c + nhpb= nhpb0 + link_start=1 + link_end=nhpb + wstrain=wstrain0 + + return + end +c-------------------------------------------------------- + subroutine softreg + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + 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) + time0=MPI_WTIME() + call minimize(etot,var,iretcode,nfun) + + write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun + time1=MPI_WTIME() + 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 + time0=MPI_WTIME() + call minimize(etot,var,iretcode,nfun) + write(iout,*)'SUMSL MASK DIST return code is',iretcode, + & ' eval ',nfun + ieval=nfun + + time1=MPI_WTIME() + 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 + + + time0=MPI_WTIME() + call minimize(etot,var,iretcode,nfun) + write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun + ieval=ieval+nfun + + time1=MPI_WTIME() + 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 + + time0=MPI_WTIME() + 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 + + time1=MPI_WTIME() + 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 + time0=MPI_WTIME() + call minimize(etot,var,iretcode,nfun) + write(iout,*)'------------------------------------------------' + write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun, + & '+ DIST eval',ieval + + time1=MPI_WTIME() + 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 + + + subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + 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) + integer jdata(5),isec(maxres) +c + jdata(1)=i1 + jdata(2)=i2 + jdata(3)=i3 + jdata(4)=i4 + jdata(5)=i5 + + call secondary2(.false.) + + do i=1,nres + isec(i)=0 + enddo + do j=1,nbfrag + do i=bfrag(1,j),bfrag(2,j) + isec(i)=1 + enddo + do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j)) + isec(i)=1 + enddo + enddo + do j=1,nhfrag + do i=hfrag(1,j),hfrag(2,j) + isec(i)=2 + enddo + enddo + +c +c cut strands at the ends +c + if (jdata(2)-jdata(1).gt.3) then + jdata(1)=jdata(1)+1 + jdata(2)=jdata(2)-1 + if (jdata(3).lt.jdata(4)) then + jdata(3)=jdata(3)+1 + jdata(4)=jdata(4)-1 + else + jdata(3)=jdata(3)-1 + jdata(4)=jdata(4)+1 + endif + endif + +cv call chainbuild +cv call etotal(energy(0)) +cv etot=energy(0) +cv write(iout,*) nnt,nct,etot +cv call write_pdb(ij*100,'first structure',etot) +cv write(iout,*) 'N16 test',(jdata(i),i=1,5) + +c------------------------ +c generate constrains +c + ishift=jdata(5)-2 + if(ishift.eq.0) ishift=-2 + nhpb0=nhpb + call chainbuild + do i=jdata(1),jdata(2) + isec(i)=-1 + if(jdata(4).gt.jdata(3))then + do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2 + isec(j)=-1 +cd print *,i,j,j+ishift + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=1000.0 + dhpb(nhpb)=DIST(i,j+ishift) + enddo + else + do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1 + isec(j)=-1 +cd print *,i,j,j+ishift + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=1000.0 + dhpb(nhpb)=DIST(i,j+ishift) + enddo + endif + enddo + + do i=nnt,nct-2 + do j=i+2,nct + if(isec(i).gt.0.or.isec(j).gt.0) then +cd print *,i,j + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=0.1 + dhpb(nhpb)=DIST(i,j) + endif + enddo + enddo + + call hpb_partition + + call geom_to_var(nvar,var) + maxfun0=maxfun + wstrain0=wstrain + maxfun=4000/5 + + do ico=1,5 + + wstrain=wstrain0/ico + +cv time0=MPI_WTIME() + 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=ieval+nfun +cv time1=MPI_WTIME() +cv write (iout,'(a,f6.2,f8.2,a)') +cv & ' Time for dist min.',time1-time0, +cv & nfun/(time1-time0),' eval/s' +cv call var_to_geom(nvar,var) +cv call chainbuild +cv call write_pdb(ij*100+ico,'dist cons',etot) + + enddo +c + nhpb=nhpb0 + call hpb_partition + wstrain=wstrain0 + maxfun=maxfun0 +c +cd print *,etot + wscloc0=wscloc + wscloc=10.0 + call sc_move(nnt,nct,100,100d0,nft_sc,etot) + wscloc=wscloc0 +cv call chainbuild +cv call etotal(energy(0)) +cv etot=energy(0) +cv call write_pdb(ij*100+10,'sc_move',etot) +cd call intout +cd print *,nft_sc,etot + + return + end + + subroutine beta_zip(i1,i2,ieval,ij) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + 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) + character*10 test + +cv call chainbuild +cv call etotal(energy(0)) +cv etot=energy(0) +cv write(test,'(2i5)') i1,i2 +cv call write_pdb(ij*100,test,etot) +cv write(iout,*) 'N17 test',i1,i2,etot,ij + +c +c generate constrains +c + nhpb0=nhpb + nhpb=nhpb+1 + ihpb(nhpb)=i1 + jhpb(nhpb)=i2 + forcon(nhpb)=1000.0 + dhpb(nhpb)=4.0 + + call hpb_partition + + call geom_to_var(nvar,var) + maxfun0=maxfun + wstrain0=wstrain + maxfun=1000/5 + + do ico=1,5 + wstrain=wstrain0/ico +cv time0=MPI_WTIME() + 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=ieval+nfun +cv time1=MPI_WTIME() +cv write (iout,'(a,f6.2,f8.2,a)') +cv & ' Time for dist min.',time1-time0, +cv & nfun/(time1-time0),' eval/s' +c do not comment the next line + call var_to_geom(nvar,var) +cv call chainbuild +cv call write_pdb(ij*100+ico,'dist cons',etot) + enddo + + nhpb=nhpb0 + call hpb_partition + wstrain=wstrain0 + maxfun=maxfun0 + +cv call etotal(energy(0)) +cv etot=energy(0) +cv write(iout,*) 'N17 test end',i1,i2,etot,ij + + + return + end