X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fbank.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fbank.F;h=5636ba0936060192348ea6f72b842fe490b1c2b0;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/bank.F b/source/unres/src_MD-M-newcorr/bank.F new file mode 100644 index 0000000..5636ba0 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/bank.F @@ -0,0 +1,1084 @@ +cc--------------------------------- + subroutine refresh_bank(ntrial) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.CONTROL' + character chacc + integer iaccn + double precision l_diff(mxio),denep + + do i=0,mxmv + do j=1,3 + nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j) + nstatnx(i,j)=0 + enddo + enddo + +c loop over all newly obtained conformations + do n=1,ntrial + chacc=' ' + iaccn=0 + nstatnx(movernx(n),1)=nstatnx(movernx(n),1)+1 +cccccccccccccccccccccccccccccccccccccccccccc +cjlee + if(iref.ne.0) then + if(rmsn(n).gt.rmscut.or.pncn(n).lt.pnccut) goto 100 + endif +cjlee + if(etot(n).gt.ebmax) goto 100 +c Find the conformation closest to the conformation n in the bank + difmin=9.d9 + do m=1,nbank + call get_diff12(dihang(1,1,1,n),bvar(1,1,1,m),l_diff(m)) + if(l_diff(m).lt.difmin) then + difmin=l_diff(m) + idmin=m + endif + enddo + + if(difmin.lt.cutdif) then +c n is redundant to idmin + if(etot(n).lt.bene(idmin)) then + if(etot(n).lt.bene(idmin)-0.01d0) then + ibank(idmin)=0 + jbank(idmin)=0 + endif + denep=bene(idmin)-etot(n) + call replace_bvar(idmin,n) +crc Update dij + do i1=1,nbank + if (i1.ne.idmin) then + dij(i1,idmin)=l_diff(i1) + dij(idmin,i1)=l_diff(i1) + endif + enddo + chacc='c' + iaccn=idmin + nstatnx(movernx(n),2)=nstatnx(movernx(n),2)+1 + if(idmin.eq.ibmax) call find_max + endif + else +c got new conformation + del_ene=0.0d0 + if(ebmax-ebmin.gt.del_ene) then + denep=ebmax-etot(n) + call replace_bvar(ibmax,n) +crc Update dij + do i1=1,nbank + if (i1.ne.ibmax) then + dij(i1,ibmax)=l_diff(i1) + dij(ibmax,i1)=l_diff(i1) + endif + enddo + chacc='f' + iaccn=ibmax + nstatnx(movernx(n),3)=nstatnx(movernx(n),3)+1 + ibank(ibmax)=0 + jbank(ibmax)=0 + call find_max + else + if(del_ene.lt.0.0001) then + write (iout,*) 'ERROR in refresh_bank: ' + write (iout,*) 'ebmax: ',ebmax + write (iout,*) 'ebmin: ',ebmin + write (iout,*) 'del_ene: ',del_ene +crc call mpi_abort(mpi_comm_world,ierror,ierrcode) + endif +cjp nbmax is never defined so condition below is always false +c if(nbank.lt.nbmax) then +c nbank=nbank+1 +c call replace_bvar(nbank,n) +c ibank(nbank)=0 +c jbank(nbank)=0 +c else + call replace_bvar(ibmax,n) + ibank(ibmax)=0 + jbank(ibmax)=0 + call find_max +c endif + endif + endif +cccccccccccccccccccccccccccccccccccccccccccc + 100 continue + if (iaccn.eq.0) then + if (iref.eq.0) then + write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5)') + & indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ', + & indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9) + else + write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5 + & ,a5,0pf4.1,a5,f3.0)') + & indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ', + & indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9), + & ' rms ',rmsn(n),' %NC ',pncn(n)*100 + endif + else + if (iref.eq.0) then + write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5, + & 1x,a1,i4,0pf8.1,0pf8.1)') + & indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ', + & indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9), + & chacc,iaccn,difmin,denep + else + write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5, + & 0pf4.1,a5,f3.0,1x,a1,i4,0pf8.1,0pf8.1)') + & indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ', + & indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9), + & ' rms ',rmsn(n),' %NC ',pncn(n)*100, + & chacc,iaccn,difmin,denep + endif + endif + enddo +c end of loop over all newly obtained conformations + do i=0,mxmv + if(nstatnx(i,1).ne.0) then + if (i.le.9) then + write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') + & '## N',i,' total=',nstatnx(i,1), + & ' close=',nstatnx(i,2),' far=',nstatnx(i,3), + & ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1) + else + write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') + & '##N',i,' total=',nstatnx(i,1), + & ' close=',nstatnx(i,2),' far=',nstatnx(i,3), + & ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1) + endif + else + if (i.le.9) then + write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') + & '## N',i,' total=',nstatnx(i,1), + & ' close=',nstatnx(i,2),' far=',nstatnx(i,3), + & ' %acc',0.0 + else + write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') + & '##N',i,' total=',nstatnx(i,1), + & ' close=',nstatnx(i,2),' far=',nstatnx(i,3), + & ' %acc',0.0 + endif + endif + enddo + call flush(iout) +crc Update dij +crc moved up, saves some get_diff12 calls +crc +crc do i1=1,nbank-1 +crc do i2=i1+1,nbank +crc if(jbank(i1).eq.0.or.jbank(i2).eq.0) then +crc call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff) +crc dij(i1,i2)=diff +crc dij(i2,i1)=diff +crc endif +crc enddo +crc enddo + + do i=1,nbank + jbank(i)=1 + enddo + + return + end +c--------------------------------- + subroutine replace_bvar(iold,inew) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.IOUNITS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.CHAIN' + include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' + + if (iold.gt.mxio .or. iold.lt.1 .or. inew.gt.mxio .or. inew.lt.1) + & then + write (iout,*) 'Dimension ERROR in REPLACE_BVAR: IOLD',iold, + & ' INEW',inew + call mpi_abort(mpi_comm_world,ierror,ierrcode) + endif + do k=1,numch + do j=2,nres-1 + do i=1,4 + bvar(i,j,k,iold)=dihang(i,j,k,inew) + enddo + enddo + enddo + bene(iold)=etot(inew) + brmsn(iold)=rmsn(inew) + bpncn(iold)=pncn(inew) + + if(bene(iold).lt.ebmin) then + ebmin=bene(iold) + ibmin=iold + endif + + if(vdisulf) then + bvar_nss(iold)=nss_out(inew) +cd write(iout,*) 'SS BANK',iold,bvar_nss(iold) + do i=1,bvar_nss(iold) + bvar_ss(1,i,iold)=iss_out(i,inew) + bvar_ss(2,i,iold)=jss_out(i,inew) +cd write(iout,*) 'SS',bvar_ss(1,i,iold)-nres, +cd & bvar_ss(2,i,iold)-nres + enddo + + bvar_ns(iold)=ns-2*bvar_nss(iold) +cd write(iout,*) 'CYS #free ', bvar_ns(iold) + k=0 + do i=1,ns + j=1 + do while( iss(i).ne.iss_out(j,inew)-nres .and. + & iss(i).ne.jss_out(j,inew)-nres .and. + & j.le.nss_out(inew)) + j=j+1 + enddo + if (j.gt.nss_out(inew)) then + k=k+1 + bvar_s(k,iold)=iss(i) + endif + enddo +cd write(iout,*) 'CYS free',(bvar_s(k,iold),k=1,bvar_ns(iold)) + endif + + return + end +c--------------------------------------- + subroutine write_rbank(jlee,adif,nft) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + + open(icsa_rbank,file=csa_rbank,status="unknown") + write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif + do k=1,nbank + write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k) + do j=1,numch + do l=2,nres-1 + write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4) + enddo + enddo + enddo + close(icsa_rbank) + + 850 format (10f8.3) + 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =", + & i8,i10,i2,f15.5) + 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3 + & ,' %NC ',0pf5.2) + + return + end +c--------------------------------------- + subroutine read_rbank(jlee,adif) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.IOUNITS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + include 'COMMON.SETUP' + character*80 karta + + open(icsa_rbank,file=csa_rbank,status="old") + read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif + print *,jleer,nbankr,nstepr,nftr,icycler,adif +c print *, 'adif from read_rbank ',adif + if(nbankr.ne.nbank) then + write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr, + & ' NBANK',nbank + call mpi_abort(mpi_comm_world,ierror,ierrcode) + endif + if(jleer.ne.jlee) then + write (iout,*) 'ERROR in READ_BANK: JLEER',jleer, + & ' JLEE',jlee + call mpi_abort(mpi_comm_world,ierror,ierrcode) + endif + + kk=0 + do k=1,nbankr + read (icsa_rbank,'(a80)') karta + write(iout,*) "READ_RBANK: kk=",kk + write(iout,*) karta +c if (index(karta,"*").gt.0) then +c write (iout,*) "***** Stars in bankr ***** k=",k, +c & " skipped" +c do j=1,numch +c do l=2,nres-1 +c read (30,850) (rdummy,i=1,4) +c enddo +c enddo +c else + kk=kk+1 + call reada(karta,"total E",rene(kk),1.0d20) + call reada(karta,"rmsd from N",rrmsn(kk),0.0d0) + call reada(karta,"%NC",rpncn(kk),0.0d0) + write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk), + & "%NC",bpncn(kk),ibank(kk) +c read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk) + do j=1,numch + do l=2,nres-1 + read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4) +c write (iout,850) (rvar(i,l,j,kk),i=1,4) + do i=1,4 + rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk) + enddo + enddo + enddo +c endif + enddo +cd write (*,*) "read_rbank ******************* kk",kk, +cd & "nbankr",nbankr + if (kk.lt.nbankr) nbankr=kk +cd do kk=1,nbankr +cd print *,"kk=",kk +cd do j=1,numch +cd do l=2,nres-1 +cd write (*,850) (rvar(i,l,j,kk),i=1,4) +cd enddo +cd enddo +cd enddo + close(icsa_rbank) + + 850 format (10f8.3) + 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5) + 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2) + + return + end +c--------------------------------------- + subroutine write_bank(jlee,nft) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + include 'COMMON.SBRIDGE' + include 'COMMON.CONTROL' + character*7 chtmp + character*40 chfrm + external ilen + + open(icsa_bank,file=csa_bank,status="unknown") + write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif + write (icsa_bank,902) nglob_csa, eglob_csa + open (igeom,file=intname,status='UNKNOWN') + do k=1,nbank + write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k) + if (vdisulf) write (icsa_bank,'(101i4)') + & bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k)) + do j=1,numch + do l=2,nres-1 + write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4) + enddo + enddo + if (bvar_nss(k).le.9) then + write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k), + & bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k)) + else + write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k), + & bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9) + write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k), + & bvar_ss(2,i,k),i=10,bvar_nss(k)) + endif + write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1) + write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2) + write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1) + write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1) + enddo + close(icsa_bank) + close(igeom) + + if (nstep/200.gt.ilastnstep) then + + ilastnstep=(ilastnstep+1)*1.5 + write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')' + write(chtmp,chfrm) nstep + open(icsa_int,file=prefix(:ilen(prefix)) + & //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN') + do k=1,nbank + if (bvar_nss(k).le.9) then + write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k), + & bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k)) + else + write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k), + & bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9) + write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k), + & bvar_ss(2,i,k),i=10,bvar_nss(k)) + endif + write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1) + write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2) + write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1) + write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1) + enddo + close(icsa_int) + endif + + + 200 format (8f10.4) + 850 format (10f8.3) + 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =", + & i8,i10,i2,f15.5) + 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5) + 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3, + & ' %NC ',0pf5.2,i5) + + return + end +c--------------------------------------- + subroutine write_bank_reminimized(jlee,nft) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + include 'COMMON.SBRIDGE' + + open(icsa_bank_reminimized,file=csa_bank_reminimized, + & status="unknown") + write (icsa_bank_reminimized,900) + & jlee,nbank,nstep,nft,icycle,cutdif + open (igeom,file=intname,status='UNKNOWN') + do k=1,nbank + write (icsa_bank_reminimized,952) k,bene(k),brmsn(k), + & bpncn(k),ibank(k) + do j=1,numch + do l=2,nres-1 + write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4) + enddo + enddo + if (nss.le.9) then + write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k), + & nss,(ihpb(i),jhpb(i),i=1,nss) + else + write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k), + & nss,(ihpb(i),jhpb(i),i=1,9) + write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss) + endif + write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1) + write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2) + write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1) + write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1) + enddo + close(icsa_bank_reminimized) + close(igeom) + + 200 format (8f10.4) + 850 format (10f8.3) + 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =", + & i8,i10,i2,f15.5) + 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3 + & ,' %NC ',0pf5.2,i5) + + return + end +c--------------------------------- + subroutine read_bank(jlee,nft,cutdifr) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' + character*80 karta + integer ilen + external ilen + + open(icsa_bank,file=csa_bank,status="old") + read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr + read (icsa_bank,902) nglob_csa, eglob_csa +c if(jleer.ne.jlee) then +c write (iout,*) 'ERROR in READ_BANK: JLEER',jleer, +c & ' JLEE',jlee +c call mpi_abort(mpi_comm_world,ierror,ierrcode) +c endif + + kk=0 + do k=1,nbank + read (icsa_bank,'(a80)') karta + write(iout,*) "READ_BANK: kk=",kk + write(iout,*) karta +c if (index(karta,"*").gt.0) then +c write (iout,*) "***** Stars in bank ***** k=",k, +c & " skipped" +c do j=1,numch +c do l=2,nres-1 +c read (33,850) (rdummy,i=1,4) +c enddo +c enddo +c else + kk=kk+1 + call reada(karta,"total E",bene(kk),1.0d20) + call reada(karta,"rmsd from N",brmsn(kk),0.0d0) + call reada(karta,"%NC",bpncn(kk),0.0d0) + read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk) + goto 112 + 111 ibank(kk)=0 + 112 continue + write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk), + & "%NC",bpncn(kk),ibank(kk) +c read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k) + if (vdisulf) then + read (icsa_bank,'(101i4)') + & bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk)) + bvar_ns(kk)=ns-2*bvar_nss(kk) + write(iout,*) 'read SSBOND',bvar_nss(kk), + & ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk)) +cd write(iout,*) 'read CYS #free ', bvar_ns(kk) + l=0 + do i=1,ns + j=1 + do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. + & iss(i).ne.bvar_ss(2,j,kk)-nres .and. + & j.le.bvar_nss(kk)) + j=j+1 + enddo + if (j.gt.bvar_nss(kk)) then + l=l+1 + bvar_s(l,kk)=iss(i) + endif + enddo +cd write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk)) + endif + do j=1,numch + do l=2,nres-1 + read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4) +c write (iout,850) (bvar(i,l,j,kk),i=1,4) + do i=1,4 + bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk) + enddo ! l + enddo ! l + enddo ! j +c endif + enddo ! k + + if (kk.lt.nbank) nbank=kk +cd write (*,*) "read_bank ******************* kk",kk, +cd & "nbank",nbank +cd do kk=1,nbank +cd print *,"kk=",kk +cd do j=1,numch +cd do l=2,nres-1 +cd write (*,850) (bvar(i,l,j,kk),i=1,4) +cd enddo +cd enddo +cd enddo + +c do k=1,nbank +c read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k) +c do j=1,numch +c do l=2,nres-1 +c read (33,850) (bvar(i,l,j,k),i=1,4) +c do i=1,4 +c bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k) +c enddo +c enddo +c enddo +c enddo + close(icsa_bank) + + 850 format (10f8.3) + 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5) + 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5) + 902 format (1x,11x,i4,12x,1pe14.5) + 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5) + + return + end +c--------------------------------------- + subroutine write_bank1(jlee) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.CHAIN' + include 'COMMON.GEO' + +#if defined(AIX) || defined(PGI) + open(icsa_bank1,file=csa_bank1,position="append") +#else + open(icsa_bank1,file=csa_bank1,access="append") +#endif + write (icsa_bank1,900) jlee,nbank,nstep,cutdif + do k=1,nbank + write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k) + do j=1,numch + do l=2,nres-1 + write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4) + enddo + enddo + enddo + close(icsa_bank1) + 850 format (10f8.3) + 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5) + 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3 + & ,' %NC ',0pf5.2,i5) + + return + end +c--------------------------------- + subroutine save_is(ind) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.IOUNITS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.CHAIN' + + index=nbank+ind +c print *, "nbank,ind,index,is(ind) ",nbank,ind,index,is(ind) + if (index.gt.mxio .or. index.lt.1 .or. + & is(ind).gt.mxio .or. is(ind).lt.1) then + write (iout,*) 'Dimension ERROR in SAVE_IS: INDEX',index, + & ' IND',ind,' IS',is(ind) + call mpi_abort(mpi_comm_world,ierror,ierrcode) + endif + do k=1,numch + do j=2,nres-1 + do i=1,4 + bvar(i,j,k,index)=bvar(i,j,k,is(ind)) + enddo + enddo + enddo + bene(index)=bene(is(ind)) + ibank(is(ind))=1 + + return + end +c--------------------------------- + subroutine select_is(n,ifar,idum) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + dimension itag(mxio),adiff(mxio) + + iuse=0 + do i=1,nbank + if(ibank(i).eq.0) then + iuse=iuse+1 + itag(iuse)=i + endif + enddo + iusesv=iuse + + if(iuse.eq.0) then + icycle=icycle+1 + do i=1,nbank + if(ibank(i).eq.2) then + ibank(i)=1 + else + ibank(i)=0 + endif + enddo + imade=0 + call get_is(idum,ifar,n,imade,0) +ctest3 call get_is_max(idum,ifar,n,imade,0) + else if(iuse.eq.n) then + do i=1,iuse + is(i)=itag(i) + call save_is(i) + enddo + else if(iuse.lt.n) then +c if(icycle.eq.0) then +c do i=1,n +c ind=mod(i-1,iuse)+1 +c is(i)=itag(ind) +c call save_is(i) +c enddo +c else +c endif + do i=1,iuse + is(i)=itag(i) + call save_is(i) + enddo + imade=iuse +c call get_is_ran(idum,n,imade,1) + call get_is(idum,ifar,n,imade,1) +ctest3 call get_is_max(idum,ifar,n,imade,1) +c if(iusesv.le.n/10) then + if(iusesv.le.0) then + icycle=icycle+1 + do i=1,nbank +c if(ibank(i).eq.2) then +c ibank(i)=1 + if(ibank(i).ge.2) then + ibank(i)=ibank(i)-1 + else + ibank(i)=0 + endif + enddo + endif + else + imade=0 + call get_is(idum,ifar,n,imade,0) +ctest3 call get_is_max(idum,ifar,n,imade,0) + endif + iuse=iusesv + + return + end +c--------------------------------- + subroutine get_is_ran(idum,n,imade,k) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + real ran1,ran2 + dimension itag(mxio),adiff(mxio) + + do j=imade+1,n + iuse=0 + do i=1,nbank + if(ibank(i).eq.k) then + iuse=iuse+1 + itag(iuse)=i + endif + enddo + iran=iuse* ran1(idum)+1 + is(j)=itag(iran) + call save_is(j) + enddo + + return + end +c--------------------------------- + subroutine get_is(idum,ifar,n,imade,k) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + real ran1,ran2 + dimension itag(mxio),adiff(mxio) + + iuse=0 + do i=1,nbank + if(ibank(i).eq.k) then + iuse=iuse+1 + itag(iuse)=i + endif + enddo + iran=iuse* ran1(idum)+1 + imade=imade+1 + is(imade)=itag(iran) + call save_is(imade) + + do i=imade+1,ifar-1 + if(icycle.eq.-1) then + call select_iseed_max(i,k) + else + call select_iseed_min(i,k) +ctest4 call select_iseed_max(i,k) + endif + call save_is(i) + enddo + + do i=ifar,n + call select_iseed_far(i,k) + call save_is(i) + enddo + + return + end +c--------------------------------- + subroutine select_iseed_max(imade1,ik) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + dimension itag(mxio),adiff(mxio) + + iuse=0 + avedif=0.d0 + difmax=0.d0 + do n=1,nbank + if(ibank(n).eq.ik) then + iuse=iuse+1 + diffmn=9.d190 + do imade=1,imade1-1 +c m=nbank+imade +c call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff) + m=is(imade) + diff=dij(n,m) + if(diff.lt.diffmn) diffmn=diff + enddo + if(diffmn.gt.difmax) difmax=diffmn + adiff(iuse)=diffmn + itag(iuse)=n + avedif=avedif+diffmn + endif + enddo + + avedif=avedif/iuse +c avedif=(avedif+difmax)/2 + emax=-9.d190 + do i=1,iuse + if(adiff(i).ge.avedif) then + itagi=itag(i) + benei=bene(itagi) + if(benei.gt.emax) then + emax=benei + is(imade1)=itagi + endif + endif + enddo + + if(ik.eq.0) iuse=iuse-1 + + return + end +c--------------------------------- + subroutine select_iseed_min(imade1,ik) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + dimension itag(mxio),adiff(mxio) + + iuse=0 + avedif=0.d0 + difmax=0.d0 + do n=1,nbank + if(ibank(n).eq.ik) then + iuse=iuse+1 + diffmn=9.d190 + do imade=1,imade1-1 +c m=nbank+imade +c call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff) + m=is(imade) + diff=dij(n,m) + if(diff.lt.diffmn) diffmn=diff + enddo + if(diffmn.gt.difmax) difmax=diffmn + adiff(iuse)=diffmn + itag(iuse)=n + avedif=avedif+diffmn + endif + enddo + + avedif=avedif/iuse +c avedif=(avedif+difmax)/2 + emin=9.d190 + do i=1,iuse +c print *,"i, adiff(i),avedif : ",i,adiff(i),avedif + if(adiff(i).ge.avedif) then + itagi=itag(i) + benei=bene(itagi) +c print *,"i, benei,emin : ",i,benei,emin + if(benei.lt.emin) then + emin=benei + is(imade1)=itagi + endif + endif + enddo + + if(ik.eq.0) iuse=iuse-1 + +c print *, "exiting select_iseed_min",is(imade1) + + return + end +c--------------------------------- + subroutine select_iseed_far(imade1,ik) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + + dmax=-9.d190 + do n=1,nbank + if(ibank(n).eq.ik) then + diffmn=9.d190 + do imade=1,imade1-1 +c m=nbank+imade +c call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff) + m=is(imade) + diff=dij(n,m) + if(diff.lt.diffmn) diffmn=diff + enddo + endif + if(diffmn.gt.dmax) then + dmax=diffmn + is(imade1)=n + endif + enddo + + return + end +c--------------------------------- + subroutine find_min + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + + ebmin=9.d190 + + do i=1,nbank + benei=bene(i) + if(benei.lt.ebmin) then + ebmin=benei + ibmin=i + endif + enddo + + return + end +c--------------------------------- + subroutine write_csa_pdb(var,ene,nft,ik,iw_pdb) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.MINIM' + include 'COMMON.SETUP' + include 'COMMON.GEO' + include 'COMMON.CHAIN' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.SBRIDGE' + integer lenpre,lenpot,ilen + external ilen + dimension var(maxvar) + character*50 titelloc + character*3 zahl + + nmin_csa=nmin_csa+1 + if(ene.lt.eglob_csa) then + eglob_csa=ene + nglob_csa=nglob_csa+1 + call numstr(nglob_csa,zahl) + + call var_to_geom(nvar,var) + call chainbuild + call secondary2(.false.) + + lenpre=ilen(prefix) + open(icsa_pdb,file=prefix(:lenpre)//'@'//zahl//'.pdb') + + if (iw_pdb.eq.1) then + write(titelloc,'(a2,i3,a3,i9,a3,i6)') + & 'GM',nglob_csa,' e ',nft,' m ',nmin_csa + else + write(titelloc,'(a2,i3,a3,i9,a3,i6,a5,f5.2,a5,f5.1)') + & 'GM',nglob_csa,' e ',nft,' m ',nmin_csa,' rms ' + & ,rmsn(ik),' %NC ',pncn(ik)*100 + endif + call pdbout(eglob_csa,titelloc,icsa_pdb) + close(icsa_pdb) + endif + + return + end +c--------------------------------- + subroutine find_max + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + + ebmax=-9.d190 + + do i=1,nbank + benei=bene(i) + if(benei.gt.ebmax) then + ebmax=benei + ibmax=i + endif + enddo + + return + end +c--------------------------------- + subroutine get_diff + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + + tdiff=0.d0 + difmin=9.d190 + do i1=1,nbank-1 + do i2=i1+1,nbank + if(jbank(i1).eq.0.or.jbank(i2).eq.0) then + call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff) + dij(i1,i2)=diff + dij(i2,i1)=diff + else + diff=dij(i1,i2) + endif + tdiff=tdiff+diff + if(diff.lt.difmin) difmin=diff + enddo + dij(i1,i1)=0.0 + enddo + + do i=1,nbank + jbank(i)=1 + enddo + + avedif=tdiff/nbank/(nbank-1)*2 + + return + end +c--------------------------------- + subroutine estimate_cutdif(adif,xct,cutdifr) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + + ctdif1=adif/cut2 + + exponent = cutdifr*cut1/adif + exponent = dlog(exponent)/dlog(xct) + + nexp=exponent+0.25 + cutdif= adif/cut1*xct**nexp + if(cutdif.lt.ctdif1) cutdif=ctdif1 + + return + end +c--------------------------------- + subroutine get_is_max(idum,ifar,n,imade,k) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CSA' + include 'COMMON.BANK' + double precision emax + + do i=imade+1,n + emax=-9.d190 + do j=1,nbank + if(ibank(j).eq.k .and. bene(j).gt.emax) then + emax=bene(j) + is(i)=j + endif + enddo + call save_is(i) + enddo + + return + end