source/unres/src_CSA/bank.F

   1 cc---------------------------------
   2       subroutine refresh_bank(ntrial)
   3       implicit real*8 (a-h,o-z)
   4       include 'DIMENSIONS'
   5       include 'mpif.h'
   6       include 'COMMON.CSA'
   7       include 'COMMON.BANK'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.CHAIN'
  10       include 'COMMON.VAR'
  11       include 'COMMON.CONTROL'
  12       character chacc
  13       integer iaccn
  14       double precision l_diff(mxio),denep
  15
  16       do i=0,mxmv
  17         do j=1,3
  18           nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j)
  19           nstatnx(i,j)=0
  20         enddo
  21       enddo
  22
  23 c loop over all newly obtained conformations
  24       do n=1,ntrial
  25        chacc=' '
  26        iaccn=0
  27        nstatnx(movernx(n),1)=nstatnx(movernx(n),1)+1
  28 cccccccccccccccccccccccccccccccccccccccccccc
  29 cjlee
  30        if(iref.ne.0) then
  31         if(rmsn(n).gt.rmscut.or.pncn(n).lt.pnccut) goto 100
  32        endif
  33 cjlee
  34        if(etot(n).gt.ebmax) goto 100
  35 c Find the conformation closest to the conformation n in the bank
  36        difmin=9.d9
  37        do m=1,nbank
  38         call get_diff12(dihang(1,1,1,n),bvar(1,1,1,m),l_diff(m))
  39         if(l_diff(m).lt.difmin) then
  40          difmin=l_diff(m)
  41          idmin=m
  42         endif
  43        enddo
  44
  45        if(difmin.lt.cutdif) then
  46 c n is redundant to idmin
  47         if(etot(n).lt.bene(idmin)) then
  48          if(etot(n).lt.bene(idmin)-0.01d0) then
  49           ibank(idmin)=0
  50           jbank(idmin)=0
  51          endif
  52          denep=bene(idmin)-etot(n)
  53          call replace_bvar(idmin,n)
  54 crc Update dij
  55          do i1=1,nbank
  56            if (i1.ne.idmin) then
  57             dij(i1,idmin)=l_diff(i1)
  58             dij(idmin,i1)=l_diff(i1)
  59            endif
  60          enddo
  61          chacc='c'
  62          iaccn=idmin
  63          nstatnx(movernx(n),2)=nstatnx(movernx(n),2)+1
  64          if(idmin.eq.ibmax) call find_max
  65         endif
  66        else
  67 c got new conformation
  68         del_ene=0.0d0
  69         if(ebmax-ebmin.gt.del_ene) then
  70          denep=ebmax-etot(n)
  71          call replace_bvar(ibmax,n)
  72 crc Update dij
  73          do i1=1,nbank
  74            if (i1.ne.ibmax) then
  75             dij(i1,ibmax)=l_diff(i1)
  76             dij(ibmax,i1)=l_diff(i1)
  77            endif
  78          enddo
  79          chacc='f'
  80          iaccn=ibmax
  81          nstatnx(movernx(n),3)=nstatnx(movernx(n),3)+1
  82          ibank(ibmax)=0
  83          jbank(ibmax)=0
  84          call find_max
  85         else
  86          if(del_ene.lt.0.0001) then
  87           write (iout,*) 'ERROR in refresh_bank: '
  88           write (iout,*) 'ebmax: ',ebmax
  89           write (iout,*) 'ebmin: ',ebmin
  90           write (iout,*) 'del_ene: ',del_ene
  91 crc          call mpi_abort(mpi_comm_world,ierror,ierrcode)
  92          endif
  93 cjp nbmax is never defined so condition below is always false
  94 c         if(nbank.lt.nbmax) then
  95 c          nbank=nbank+1
  96 c          call replace_bvar(nbank,n)
  97 c          ibank(nbank)=0
  98 c          jbank(nbank)=0
  99 c         else
 100           call replace_bvar(ibmax,n)
 101           ibank(ibmax)=0
 102           jbank(ibmax)=0
 103           call find_max
 104 c         endif
 105         endif
 106        endif
 107 cccccccccccccccccccccccccccccccccccccccccccc
 108   100 continue
 109        if (iaccn.eq.0) then
 110         if (iref.eq.0) then
 111          write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5)')
 112      &    indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',
 113      &    indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9)
 114         else
 115       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5
 116      &       ,a5,0pf4.1,a5,f3.0)')
 117      &    indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',
 118      &    indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),
 119      &    ' rms ',rmsn(n),' %NC ',pncn(n)*100
 120         endif
 121        else
 122         if (iref.eq.0) then
 123         write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,
 124      &         1x,a1,i4,0pf8.1,0pf8.1)')
 125      &    indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',
 126      &    indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),
 127      &    chacc,iaccn,difmin,denep
 128         else
 129          write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5,
 130      &                   0pf4.1,a5,f3.0,1x,a1,i4,0pf8.1,0pf8.1)')
 131      &    indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',
 132      &    indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),
 133      &    ' rms ',rmsn(n),' %NC ',pncn(n)*100,
 134      &     chacc,iaccn,difmin,denep
 135         endif
 136        endif
 137       enddo
 138 c end of loop over all newly obtained conformations
 139       do i=0,mxmv
 140         if(nstatnx(i,1).ne.0) then
 141          if (i.le.9) then
 142          write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)')
 143      &          '## N',i,' total=',nstatnx(i,1),
 144      &          ' close=',nstatnx(i,2),' far=',nstatnx(i,3),
 145      &          ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1)
 146          else
 147          write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)')
 148      &          '##N',i,' total=',nstatnx(i,1),
 149      &          ' close=',nstatnx(i,2),' far=',nstatnx(i,3),
 150      &          ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1)
 151          endif
 152         else
 153          if (i.le.9) then
 154          write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)')
 155      &          '## N',i,' total=',nstatnx(i,1),
 156      &          ' close=',nstatnx(i,2),' far=',nstatnx(i,3),
 157      &          ' %acc',0.0
 158          else
 159          write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)')
 160      &          '##N',i,' total=',nstatnx(i,1),
 161      &          ' close=',nstatnx(i,2),' far=',nstatnx(i,3),
 162      &          ' %acc',0.0
 163          endif
 164         endif
 165       enddo
 166       call flush(iout)
 167 crc Update dij
 168 crc moved up, saves some get_diff12 calls
 169 crc
 170 crc      do i1=1,nbank-1
 171 crc       do i2=i1+1,nbank
 172 crc        if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
 173 crc         call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
 174 crc         dij(i1,i2)=diff
 175 crc         dij(i2,i1)=diff
 176 crc        endif
 177 crc       enddo
 178 crc      enddo
 179
 180       do i=1,nbank
 181        jbank(i)=1
 182       enddo
 183
 184       return
 185       end
 186 c---------------------------------
 187       subroutine replace_bvar(iold,inew)
 188       implicit real*8 (a-h,o-z)
 189       include 'DIMENSIONS'
 190       include 'mpif.h'
 191       include 'COMMON.IOUNITS'
 192       include 'COMMON.CSA'
 193       include 'COMMON.BANK'
 194       include 'COMMON.CHAIN'
 195       include 'COMMON.CONTROL'
 196       include 'COMMON.SBRIDGE'
 197
 198       if (iold.gt.mxio .or. iold.lt.1 .or. inew.gt.mxio .or. inew.lt.1)
 199      &  then
 200         write (iout,*) 'Dimension ERROR in REPLACE_BVAR: IOLD',iold,
 201      &  ' INEW',inew
 202         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 203       endif
 204       do k=1,numch
 205        do j=2,nres-1
 206         do i=1,4
 207          bvar(i,j,k,iold)=dihang(i,j,k,inew)
 208         enddo
 209        enddo
 210       enddo
 211       bene(iold)=etot(inew)
 212       brmsn(iold)=rmsn(inew)
 213       bpncn(iold)=pncn(inew)
 214
 215       if(bene(iold).lt.ebmin) then
 216         ebmin=bene(iold)
 217         ibmin=iold
 218       endif
 219
 220       if(vdisulf) then
 221         bvar_nss(iold)=nss_out(inew)
 222 cd        write(iout,*) 'SS BANK',iold,bvar_nss(iold)
 223         do i=1,bvar_nss(iold)
 224           bvar_ss(1,i,iold)=iss_out(i,inew)
 225           bvar_ss(2,i,iold)=jss_out(i,inew)
 226 cd          write(iout,*) 'SS',bvar_ss(1,i,iold)-nres,
 227 cd     &          bvar_ss(2,i,iold)-nres
 228         enddo
 229
 230         bvar_ns(iold)=ns-2*bvar_nss(iold)
 231 cd        write(iout,*) 'CYS #free ', bvar_ns(iold)
 232            k=0
 233            do i=1,ns
 234              j=1
 235              do while( iss(i).ne.iss_out(j,inew)-nres .and.
 236      &                 iss(i).ne.jss_out(j,inew)-nres .and.
 237      &                 j.le.nss_out(inew))
 238                 j=j+1
 239              enddo
 240              if (j.gt.nss_out(inew)) then
 241                k=k+1
 242                bvar_s(k,iold)=iss(i)
 243              endif
 244            enddo
 245 cd         write(iout,*) 'CYS free',(bvar_s(k,iold),k=1,bvar_ns(iold))
 246       endif
 247
 248       return
 249       end
 250 c---------------------------------------
 251       subroutine write_rbank(jlee,adif,nft)
 252       implicit real*8 (a-h,o-z)
 253       include 'DIMENSIONS'
 254       include 'COMMON.IOUNITS'
 255       include 'COMMON.CSA'
 256       include 'COMMON.BANK'
 257       include 'COMMON.CHAIN'
 258       include 'COMMON.GEO'
 259
 260       open(icsa_rbank,file=csa_rbank,status="unknown")
 261       write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
 262       do k=1,nbank
 263        write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
 264        do j=1,numch
 265         do l=2,nres-1
 266          write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
 267         enddo
 268        enddo
 269       enddo
 270       close(icsa_rbank)
 271
 272   850 format (10f8.3)
 273   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",
 274      &        i8,i10,i2,f15.5)
 275   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3
 276      &           ,' %NC ',0pf5.2)
 277
 278       return
 279       end
 280 c---------------------------------------
 281       subroutine read_rbank(jlee,adif)
 282       implicit real*8 (a-h,o-z)
 283       include 'DIMENSIONS'
 284       include 'mpif.h'
 285       include 'COMMON.IOUNITS'
 286       include 'COMMON.CSA'
 287       include 'COMMON.BANK'
 288       include 'COMMON.CHAIN'
 289       include 'COMMON.GEO'
 290       include 'COMMON.SETUP'
 291       character*80 karta
 292
 293       open(icsa_rbank,file=csa_rbank,status="old")
 294       read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
 295       print *,jleer,nbankr,nstepr,nftr,icycler,adif
 296 c       print *, 'adif from read_rbank ',adif
 297       if(nbankr.ne.nbank) then
 298         write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,
 299      &  ' NBANK',nbank
 300         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 301       endif
 302       if(jleer.ne.jlee) then
 303         write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
 304      &  ' JLEE',jlee
 305         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 306       endif
 307
 308       kk=0
 309       do k=1,nbankr
 310         read (icsa_rbank,'(a80)') karta
 311         write(iout,*) "READ_RBANK: kk=",kk
 312         write(iout,*) karta
 313 c        if (index(karta,"*").gt.0) then
 314 c          write (iout,*) "***** Stars in bankr ***** k=",k,
 315 c     &      " skipped"
 316 c          do j=1,numch
 317 c            do l=2,nres-1
 318 c              read (30,850) (rdummy,i=1,4)
 319 c            enddo
 320 c          enddo
 321 c        else
 322           kk=kk+1
 323           call reada(karta,"total E",rene(kk),1.0d20)
 324           call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
 325           call reada(karta,"%NC",rpncn(kk),0.0d0)
 326           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),
 327      &      "%NC",bpncn(kk),ibank(kk)
 328 c          read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
 329           do j=1,numch
 330             do l=2,nres-1
 331               read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
 332 c              write (iout,850) (rvar(i,l,j,kk),i=1,4)
 333               do i=1,4
 334                 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
 335               enddo
 336             enddo
 337           enddo
 338 c        endif
 339       enddo
 340 cd      write (*,*) "read_rbank ******************* kk",kk,
 341 cd     &  "nbankr",nbankr
 342       if (kk.lt.nbankr) nbankr=kk
 343 cd      do kk=1,nbankr
 344 cd          print *,"kk=",kk
 345 cd          do j=1,numch
 346 cd            do l=2,nres-1
 347 cd              write (*,850) (rvar(i,l,j,kk),i=1,4)
 348 cd            enddo
 349 cd          enddo
 350 cd      enddo
 351       close(icsa_rbank)
 352
 353   850 format (10f8.3)
 354   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 355   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
 356
 357       return
 358       end
 359 c---------------------------------------
 360       subroutine write_bank(jlee,nft)
 361       implicit real*8 (a-h,o-z)
 362       include 'DIMENSIONS'
 363       include 'COMMON.IOUNITS'
 364       include 'COMMON.CSA'
 365       include 'COMMON.BANK'
 366       include 'COMMON.CHAIN'
 367       include 'COMMON.GEO'
 368       include 'COMMON.SBRIDGE'
 369       include 'COMMON.CONTROL'
 370       character*7 chtmp
 371       character*40 chfrm
 372       external ilen
 373
 374       open(icsa_bank,file=csa_bank,status="unknown")
 375       write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
 376       write (icsa_bank,902) nglob_csa, eglob_csa
 377       open (igeom,file=intname,status='UNKNOWN')
 378       do k=1,nbank
 379        write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 380        if (vdisulf) write (icsa_bank,'(101i4)')
 381      &    bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
 382        do j=1,numch
 383         do l=2,nres-1
 384          write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 385         enddo
 386        enddo
 387        if (bvar_nss(k).le.9) then
 388          write (igeom,'(I5,F10.2,I2,9(1X,2I3))') k,bene(k),
 389      &     bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 390        else
 391          write (igeom,'(I5,F10.2,I2,9(1X,2I3))') k,bene(k),
 392      &     bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 393          write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),
 394      &                                bvar_ss(2,i,k),i=10,bvar_nss(k))
 395        endif
 396        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 397        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 398        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 399        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 400       enddo
 401       close(icsa_bank)
 402       close(igeom)
 403
 404       if (nstep/200.gt.ilastnstep) then
 405
 406        ilastnstep=(ilastnstep+1)*1.5
 407        write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
 408        write(chtmp,chfrm) nstep
 409        open(icsa_int,file=prefix(:ilen(prefix))
 410      &         //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
 411        do k=1,nbank
 412         if (bvar_nss(k).le.9) then
 413          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),
 414      &  bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 415         else
 416          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),
 417      &     bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 418          write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),
 419      &                          bvar_ss(2,i,k),i=10,bvar_nss(k))
 420         endif
 421         write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 422         write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 423         write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 424         write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 425        enddo
 426        close(icsa_int)
 427       endif
 428
 429
 430   200 format (8f10.4)
 431   850 format (10f8.3)
 432   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",
 433      &        i8,i10,i2,f15.5)
 434   902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
 435   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,
 436      &        ' %NC ',0pf5.2,i5)
 437
 438       return
 439       end
 440 c---------------------------------------
 441       subroutine write_bank_reminimized(jlee,nft)
 442       implicit real*8 (a-h,o-z)
 443       include 'DIMENSIONS'
 444       include 'COMMON.IOUNITS'
 445       include 'COMMON.CSA'
 446       include 'COMMON.BANK'
 447       include 'COMMON.CHAIN'
 448       include 'COMMON.GEO'
 449       include 'COMMON.SBRIDGE'
 450
 451       open(icsa_bank_reminimized,file=csa_bank_reminimized,
 452      &  status="unknown")
 453       write (icsa_bank_reminimized,900)
 454      &  jlee,nbank,nstep,nft,icycle,cutdif
 455       open (igeom,file=intname,status='UNKNOWN')
 456       do k=1,nbank
 457        write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),
 458      &  bpncn(k),ibank(k)
 459        do j=1,numch
 460         do l=2,nres-1
 461          write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 462         enddo
 463        enddo
 464        if (nss.le.9) then
 465          write (igeom,'(I5,F10.2,I2,9(1X,2I3))') k,bene(k),
 466      &     nss,(ihpb(i),jhpb(i),i=1,nss)
 467        else
 468          write (igeom,'(I5,F10.2,I2,9(1X,2I3))') k,bene(k),
 469      &     nss,(ihpb(i),jhpb(i),i=1,9)
 470          write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
 471        endif
 472        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 473        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 474        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 475        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 476       enddo
 477       close(icsa_bank_reminimized)
 478       close(igeom)
 479
 480   200 format (8f10.4)
 481   850 format (10f8.3)
 482   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",
 483      &        i8,i10,i2,f15.5)
 484   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3
 485      &        ,' %NC ',0pf5.2,i5)
 486
 487       return
 488       end
 489 c---------------------------------
 490       subroutine read_bank(jlee,nft,cutdifr)
 491       implicit real*8 (a-h,o-z)
 492       include 'DIMENSIONS'
 493       include 'COMMON.IOUNITS'
 494       include 'COMMON.CSA'
 495       include 'COMMON.BANK'
 496       include 'COMMON.CHAIN'
 497       include 'COMMON.GEO'
 498       include 'COMMON.CONTROL'
 499       include 'COMMON.SBRIDGE'
 500       character*80 karta
 501       integer ilen
 502       external ilen
 503
 504       open(icsa_bank,file=csa_bank,status="old")
 505        read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
 506        read (icsa_bank,902) nglob_csa, eglob_csa
 507 c      if(jleer.ne.jlee) then
 508 c        write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
 509 c    &   ' JLEE',jlee
 510 c        call mpi_abort(mpi_comm_world,ierror,ierrcode)
 511 c      endif
 512
 513       kk=0
 514       do k=1,nbank
 515         read (icsa_bank,'(a80)') karta
 516         write(iout,*) "READ_BANK: kk=",kk
 517         write(iout,*) karta
 518 c        if (index(karta,"*").gt.0) then
 519 c          write (iout,*) "***** Stars in bank ***** k=",k,
 520 c     &      " skipped"
 521 c          do j=1,numch
 522 c            do l=2,nres-1
 523 c              read (33,850) (rdummy,i=1,4)
 524 c            enddo
 525 c          enddo
 526 c        else
 527           kk=kk+1
 528           call reada(karta,"total E",bene(kk),1.0d20)
 529           call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
 530           call reada(karta,"%NC",bpncn(kk),0.0d0)
 531           read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
 532           goto 112
 533   111     ibank(kk)=0
 534   112     continue
 535           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),
 536      &      "%NC",bpncn(kk),ibank(kk)
 537 c          read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 538           if (vdisulf) then
 539             read (icsa_bank,'(101i4)')
 540      &        bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 541             bvar_ns(kk)=ns-2*bvar_nss(kk)
 542             write(iout,*) 'read SSBOND',bvar_nss(kk),
 543      &                    ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 544 cd          write(iout,*) 'read CYS #free ', bvar_ns(kk)
 545             l=0
 546             do i=1,ns
 547              j=1
 548              do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and.
 549      &                 iss(i).ne.bvar_ss(2,j,kk)-nres .and.
 550      &                 j.le.bvar_nss(kk))
 551                 j=j+1
 552              enddo
 553              if (j.gt.bvar_nss(kk)) then
 554                l=l+1
 555                bvar_s(l,kk)=iss(i)
 556              endif
 557             enddo
 558 cd            write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
 559           endif
 560           do j=1,numch
 561             do l=2,nres-1
 562               read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
 563 c              write (iout,850) (bvar(i,l,j,kk),i=1,4)
 564               do i=1,4
 565                 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
 566               enddo ! l
 567             enddo ! l
 568           enddo ! j
 569 c        endif
 570       enddo ! k
 571
 572       if (kk.lt.nbank) nbank=kk
 573 cd      write (*,*) "read_bank ******************* kk",kk,
 574 cd     &  "nbank",nbank
 575 cd      do kk=1,nbank
 576 cd          print *,"kk=",kk
 577 cd          do j=1,numch
 578 cd            do l=2,nres-1
 579 cd              write (*,850) (bvar(i,l,j,kk),i=1,4)
 580 cd            enddo
 581 cd          enddo
 582 cd      enddo
 583
 584 c       do k=1,nbank
 585 c        read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 586 c        do j=1,numch
 587 c         do l=2,nres-1
 588 c          read (33,850) (bvar(i,l,j,k),i=1,4)
 589 c          do i=1,4
 590 c           bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
 591 c          enddo
 592 c         enddo
 593 c        enddo
 594 c       enddo
 595       close(icsa_bank)
 596
 597   850 format (10f8.3)
 598   952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
 599   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 600   902 format (1x,11x,i4,12x,1pe14.5)
 601   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
 602
 603       return
 604       end
 605 c---------------------------------------
 606       subroutine write_bank1(jlee)
 607       implicit real*8 (a-h,o-z)
 608       include 'DIMENSIONS'
 609       include 'COMMON.IOUNITS'
 610       include 'COMMON.CSA'
 611       include 'COMMON.BANK'
 612       include 'COMMON.CHAIN'
 613       include 'COMMON.GEO'
 614
 615 #if defined(AIX) || defined(PGI)
 616       open(icsa_bank1,file=csa_bank1,position="append")
 617 #else
 618       open(icsa_bank1,file=csa_bank1,access="append")
 619 #endif
 620       write (icsa_bank1,900) jlee,nbank,nstep,cutdif
 621       do k=1,nbank
 622        write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 623        do j=1,numch
 624         do l=2,nres-1
 625          write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 626         enddo
 627        enddo
 628       enddo
 629       close(icsa_bank1)
 630   850 format (10f8.3)
 631   900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
 632   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3
 633      &       ,' %NC ',0pf5.2,i5)
 634
 635       return
 636       end
 637 c---------------------------------
 638       subroutine save_is(ind)
 639       implicit real*8 (a-h,o-z)
 640       include 'DIMENSIONS'
 641       include 'mpif.h'
 642       include 'COMMON.IOUNITS'
 643       include 'COMMON.CSA'
 644       include 'COMMON.BANK'
 645       include 'COMMON.CHAIN'
 646
 647       index=nbank+ind
 648 c     print *, "nbank,ind,index,is(ind) ",nbank,ind,index,is(ind)
 649       if (index.gt.mxio .or. index.lt.1 .or.
 650      &  is(ind).gt.mxio .or. is(ind).lt.1) then
 651         write (iout,*) 'Dimension ERROR in SAVE_IS: INDEX',index,
 652      &  ' IND',ind,' IS',is(ind)
 653         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 654       endif
 655       do k=1,numch
 656        do j=2,nres-1
 657         do i=1,4
 658          bvar(i,j,k,index)=bvar(i,j,k,is(ind))
 659         enddo
 660        enddo
 661       enddo
 662       bene(index)=bene(is(ind))
 663       ibank(is(ind))=1
 664
 665       return
 666       end
 667 c---------------------------------
 668       subroutine select_is(n,ifar,idum)
 669       implicit real*8 (a-h,o-z)
 670       include 'DIMENSIONS'
 671       include 'COMMON.CSA'
 672       include 'COMMON.BANK'
 673       dimension itag(mxio),adiff(mxio)
 674
 675       iuse=0
 676       do i=1,nbank
 677        if(ibank(i).eq.0) then
 678         iuse=iuse+1
 679         itag(iuse)=i
 680        endif
 681       enddo
 682       iusesv=iuse
 683
 684       if(iuse.eq.0) then
 685        icycle=icycle+1
 686        do i=1,nbank
 687         if(ibank(i).eq.2) then
 688          ibank(i)=1
 689         else
 690          ibank(i)=0
 691         endif
 692        enddo
 693        imade=0
 694        call get_is(idum,ifar,n,imade,0)
 695 ctest3       call get_is_max(idum,ifar,n,imade,0)
 696       else if(iuse.eq.n) then
 697        do i=1,iuse
 698         is(i)=itag(i)
 699         call save_is(i)
 700        enddo
 701       else if(iuse.lt.n) then
 702 c      if(icycle.eq.0) then
 703 c       do i=1,n
 704 c        ind=mod(i-1,iuse)+1
 705 c        is(i)=itag(ind)
 706 c        call save_is(i)
 707 c       enddo
 708 c      else
 709 c      endif
 710        do i=1,iuse
 711         is(i)=itag(i)
 712         call save_is(i)
 713        enddo
 714        imade=iuse
 715 c      call get_is_ran(idum,n,imade,1)
 716        call get_is(idum,ifar,n,imade,1)
 717 ctest3       call get_is_max(idum,ifar,n,imade,1)
 718 c      if(iusesv.le.n/10) then
 719        if(iusesv.le.0) then
 720         icycle=icycle+1
 721         do i=1,nbank
 722 c        if(ibank(i).eq.2) then
 723 c         ibank(i)=1
 724          if(ibank(i).ge.2) then
 725           ibank(i)=ibank(i)-1
 726          else
 727           ibank(i)=0
 728          endif
 729         enddo
 730        endif
 731       else
 732        imade=0
 733        call get_is(idum,ifar,n,imade,0)
 734 ctest3       call get_is_max(idum,ifar,n,imade,0)
 735       endif
 736       iuse=iusesv
 737
 738       if (iuse.le.iucut) then
 739          icycle=icycle+1
 740          do i=1, nbank
 741             ibank(i)=0
 742          enddo
 743       endif
 744
 745
 746       return
 747       end
 748 c---------------------------------
 749       subroutine get_is_ran(idum,n,imade,k)
 750       implicit real*8 (a-h,o-z)
 751       include 'DIMENSIONS'
 752       include 'COMMON.CSA'
 753       include 'COMMON.BANK'
 754       real ran1,ran2
 755       dimension itag(mxio),adiff(mxio)
 756
 757       do j=imade+1,n
 758        iuse=0
 759        do i=1,nbank
 760         if(ibank(i).eq.k) then
 761          iuse=iuse+1
 762          itag(iuse)=i
 763         endif
 764        enddo
 765        iran=iuse*  ran1(idum)+1
 766        is(j)=itag(iran)
 767        call save_is(j)
 768       enddo
 769
 770       return
 771       end
 772 c---------------------------------
 773       subroutine get_is(idum,ifar,n,imade,k)
 774       implicit real*8 (a-h,o-z)
 775       include 'DIMENSIONS'
 776       include 'COMMON.CSA'
 777       include 'COMMON.BANK'
 778       real ran1,ran2
 779       dimension itag(mxio),adiff(mxio)
 780
 781       iuse=0
 782       do i=1,nbank
 783        if(ibank(i).eq.k) then
 784         iuse=iuse+1
 785         itag(iuse)=i
 786        endif
 787       enddo
 788       iran=iuse*  ran1(idum)+1
 789       imade=imade+1
 790       is(imade)=itag(iran)
 791       call save_is(imade)
 792
 793       do i=imade+1,ifar-1
 794        if(icycle.eq.-1) then
 795         call select_iseed_max(i,k)
 796        else
 797         call select_iseed_min(i,k)
 798 ctest4  call select_iseed_max(i,k)
 799        endif
 800        call save_is(i)
 801       enddo
 802
 803       do i=ifar,n
 804        call select_iseed_far(i,k)
 805        call save_is(i)
 806       enddo
 807
 808       return
 809       end
 810 c---------------------------------
 811       subroutine select_iseed_max(imade1,ik)
 812       implicit real*8 (a-h,o-z)
 813       include 'DIMENSIONS'
 814       include 'COMMON.CSA'
 815       include 'COMMON.BANK'
 816       dimension itag(mxio),adiff(mxio)
 817
 818       iuse=0
 819       avedif=0.d0
 820       difmax=0.d0
 821       do n=1,nbank
 822        if(ibank(n).eq.ik) then
 823         iuse=iuse+1
 824         diffmn=9.d190
 825         do imade=1,imade1-1
 826 c        m=nbank+imade
 827 c        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 828          m=is(imade)
 829          diff=dij(n,m)
 830          if(diff.lt.diffmn) diffmn=diff
 831         enddo
 832         if(diffmn.gt.difmax) difmax=diffmn
 833         adiff(iuse)=diffmn
 834         itag(iuse)=n
 835         avedif=avedif+diffmn
 836        endif
 837       enddo
 838
 839       avedif=avedif/iuse
 840 c     avedif=(avedif+difmax)/2
 841       emax=-9.d190
 842       do i=1,iuse
 843        if(adiff(i).ge.avedif) then
 844         itagi=itag(i)
 845         benei=bene(itagi)
 846         if(benei.gt.emax) then
 847          emax=benei
 848          is(imade1)=itagi
 849         endif
 850        endif
 851       enddo
 852
 853       if(ik.eq.0) iuse=iuse-1
 854
 855       return
 856       end
 857 c---------------------------------
 858       subroutine select_iseed_min(imade1,ik)
 859       implicit real*8 (a-h,o-z)
 860       include 'DIMENSIONS'
 861       include 'COMMON.CSA'
 862       include 'COMMON.BANK'
 863       dimension itag(mxio),adiff(mxio)
 864
 865       iuse=0
 866       avedif=0.d0
 867       difmax=0.d0
 868       do n=1,nbank
 869        if(ibank(n).eq.ik) then
 870         iuse=iuse+1
 871         diffmn=9.d190
 872         do imade=1,imade1-1
 873 c        m=nbank+imade
 874 c        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 875          m=is(imade)
 876          diff=dij(n,m)
 877          if(diff.lt.diffmn) diffmn=diff
 878         enddo
 879         if(diffmn.gt.difmax) difmax=diffmn
 880         adiff(iuse)=diffmn
 881         itag(iuse)=n
 882         avedif=avedif+diffmn
 883        endif
 884       enddo
 885
 886       avedif=avedif/iuse
 887 c     avedif=(avedif+difmax)/2
 888       emin=9.d190
 889       do i=1,iuse
 890 c      print *,"i, adiff(i),avedif : ",i,adiff(i),avedif
 891        if(adiff(i).ge.avedif) then
 892         itagi=itag(i)
 893         benei=bene(itagi)
 894 c       print *,"i, benei,emin : ",i,benei,emin
 895         if(benei.lt.emin) then
 896          emin=benei
 897          is(imade1)=itagi
 898         endif
 899        endif
 900       enddo
 901
 902       if(ik.eq.0) iuse=iuse-1
 903
 904 c     print *, "exiting select_iseed_min",is(imade1)
 905
 906       return
 907       end
 908 c---------------------------------
 909       subroutine select_iseed_far(imade1,ik)
 910       implicit real*8 (a-h,o-z)
 911       include 'DIMENSIONS'
 912       include 'COMMON.CSA'
 913       include 'COMMON.BANK'
 914
 915       dmax=-9.d190
 916       do n=1,nbank
 917        if(ibank(n).eq.ik) then
 918         diffmn=9.d190
 919         do imade=1,imade1-1
 920 c        m=nbank+imade
 921 c        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 922          m=is(imade)
 923          diff=dij(n,m)
 924          if(diff.lt.diffmn) diffmn=diff
 925         enddo
 926        endif
 927        if(diffmn.gt.dmax) then
 928         dmax=diffmn
 929         is(imade1)=n
 930        endif
 931       enddo
 932
 933       return
 934       end
 935 c---------------------------------
 936       subroutine find_min
 937       implicit real*8 (a-h,o-z)
 938       include 'DIMENSIONS'
 939       include 'COMMON.CSA'
 940       include 'COMMON.BANK'
 941
 942       ebmin=9.d190
 943
 944       do i=1,nbank
 945        benei=bene(i)
 946        if(benei.lt.ebmin) then
 947         ebmin=benei
 948         ibmin=i
 949        endif
 950       enddo
 951
 952       return
 953       end
 954 c---------------------------------
 955       subroutine write_csa_pdb(var,ene,nft,ik,iw_pdb)
 956       implicit real*8 (a-h,o-z)
 957       include 'DIMENSIONS'
 958       include 'COMMON.CSA'
 959       include 'COMMON.BANK'
 960       include 'COMMON.VAR'
 961       include 'COMMON.IOUNITS'
 962       include 'COMMON.MINIM'
 963       include 'COMMON.SETUP'
 964       include 'COMMON.GEO'
 965       include 'COMMON.CHAIN'
 966       include 'COMMON.LOCAL'
 967       include 'COMMON.INTERACT'
 968       include 'COMMON.NAMES'
 969       include 'COMMON.SBRIDGE'
 970       integer lenpre,lenpot,ilen
 971       external ilen
 972       dimension var(maxvar)
 973       character*50 titelloc
 974       character*3 zahl
 975
 976       nmin_csa=nmin_csa+1
 977       if(ene.lt.eglob_csa) then
 978         eglob_csa=ene
 979         nglob_csa=nglob_csa+1
 980         call numstr(nglob_csa,zahl)
 981
 982         call var_to_geom(nvar,var)
 983         call chainbuild
 984         call secondary2(.false.)
 985
 986         lenpre=ilen(prefix)
 987         open(icsa_pdb,file=prefix(:lenpre)//'@'//zahl//'.pdb')
 988
 989         if (iw_pdb.eq.1) then
 990           write(titelloc,'(a2,i3,a3,i9,a3,i6)')
 991      &    'GM',nglob_csa,' e ',nft,' m ',nmin_csa
 992         else
 993           write(titelloc,'(a2,i3,a3,i9,a3,i6,a5,f5.2,a5,f5.1)')
 994      &   'GM',nglob_csa,' e ',nft,' m ',nmin_csa,' rms '
 995      &        ,rmsn(ik),' %NC ',pncn(ik)*100
 996         endif
 997         call pdbout(eglob_csa,titelloc,icsa_pdb)
 998         close(icsa_pdb)
 999       endif
1000
1001       return
1002       end
1003 c---------------------------------
1004       subroutine find_max
1005       implicit real*8 (a-h,o-z)
1006       include 'DIMENSIONS'
1007       include 'COMMON.CSA'
1008       include 'COMMON.BANK'
1009
1010       ebmax=-9.d190
1011
1012       do i=1,nbank
1013        benei=bene(i)
1014        if(benei.gt.ebmax) then
1015         ebmax=benei
1016         ibmax=i
1017        endif
1018       enddo
1019
1020       return
1021       end
1022 c---------------------------------
1023       subroutine get_diff
1024       implicit real*8 (a-h,o-z)
1025       include 'DIMENSIONS'
1026       include 'COMMON.CSA'
1027       include 'COMMON.BANK'
1028
1029       tdiff=0.d0
1030       difmin=9.d190
1031       do i1=1,nbank-1
1032        do i2=i1+1,nbank
1033         if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
1034          call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
1035          dij(i1,i2)=diff
1036          dij(i2,i1)=diff
1037         else
1038          diff=dij(i1,i2)
1039         endif
1040         tdiff=tdiff+diff
1041         if(diff.lt.difmin) difmin=diff
1042        enddo
1043        dij(i1,i1)=0.0
1044       enddo
1045
1046       do i=1,nbank
1047        jbank(i)=1
1048       enddo
1049
1050       avedif=tdiff/nbank/(nbank-1)*2
1051
1052       return
1053       end
1054 c---------------------------------
1055       subroutine estimate_cutdif(adif,xct,cutdifr)
1056       implicit real*8 (a-h,o-z)
1057       include 'DIMENSIONS'
1058       include 'COMMON.CSA'
1059       include 'COMMON.BANK'
1060
1061       ctdif1=adif/cut2
1062
1063       exponent = cutdifr*cut1/adif
1064       exponent = dlog(exponent)/dlog(xct)
1065
1066       nexp=exponent+0.25
1067       cutdif= adif/cut1*xct**nexp
1068       if(cutdif.lt.ctdif1) cutdif=ctdif1
1069
1070       return
1071       end
1072 c---------------------------------
1073       subroutine get_is_max(idum,ifar,n,imade,k)
1074       implicit real*8 (a-h,o-z)
1075       include 'DIMENSIONS'
1076       include 'COMMON.CSA'
1077       include 'COMMON.BANK'
1078       double precision emax
1079
1080       do i=imade+1,n
1081        emax=-9.d190
1082        do j=1,nbank
1083         if(ibank(j).eq.k .and. bene(j).gt.emax) then
1084            emax=bene(j)
1085            is(i)=j
1086         endif
1087        enddo
1088        call save_is(i)
1089       enddo
1090
1091       return
1092       end