source/unres/src_MD-M/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.3,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.3,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.3,I2,9(1X,2I3))') k,bene(k),
 466      &     nss,(ihpb(i),jhpb(i),i=1,nss)
 467        else
 468          write (igeom,'(I5,F10.3,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       return
 739       end
 740 c---------------------------------
 741       subroutine get_is_ran(idum,n,imade,k)
 742       implicit real*8 (a-h,o-z)
 743       include 'DIMENSIONS'
 744       include 'COMMON.CSA'
 745       include 'COMMON.BANK'
 746       real ran1,ran2
 747       dimension itag(mxio),adiff(mxio)
 748
 749       do j=imade+1,n
 750        iuse=0
 751        do i=1,nbank
 752         if(ibank(i).eq.k) then
 753          iuse=iuse+1
 754          itag(iuse)=i
 755         endif
 756        enddo
 757        iran=iuse*  ran1(idum)+1
 758        is(j)=itag(iran)
 759        call save_is(j)
 760       enddo
 761
 762       return
 763       end
 764 c---------------------------------
 765       subroutine get_is(idum,ifar,n,imade,k)
 766       implicit real*8 (a-h,o-z)
 767       include 'DIMENSIONS'
 768       include 'COMMON.CSA'
 769       include 'COMMON.BANK'
 770       real ran1,ran2
 771       dimension itag(mxio),adiff(mxio)
 772
 773       iuse=0
 774       do i=1,nbank
 775        if(ibank(i).eq.k) then
 776         iuse=iuse+1
 777         itag(iuse)=i
 778        endif
 779       enddo
 780       iran=iuse*  ran1(idum)+1
 781       imade=imade+1
 782       is(imade)=itag(iran)
 783       call save_is(imade)
 784
 785       do i=imade+1,ifar-1
 786        if(icycle.eq.-1) then
 787         call select_iseed_max(i,k)
 788        else
 789         call select_iseed_min(i,k)
 790 ctest4  call select_iseed_max(i,k)
 791        endif
 792        call save_is(i)
 793       enddo
 794
 795       do i=ifar,n
 796        call select_iseed_far(i,k)
 797        call save_is(i)
 798       enddo
 799
 800       return
 801       end
 802 c---------------------------------
 803       subroutine select_iseed_max(imade1,ik)
 804       implicit real*8 (a-h,o-z)
 805       include 'DIMENSIONS'
 806       include 'COMMON.CSA'
 807       include 'COMMON.BANK'
 808       dimension itag(mxio),adiff(mxio)
 809
 810       iuse=0
 811       avedif=0.d0
 812       difmax=0.d0
 813       do n=1,nbank
 814        if(ibank(n).eq.ik) then
 815         iuse=iuse+1
 816         diffmn=9.d190
 817         do imade=1,imade1-1
 818 c        m=nbank+imade
 819 c        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 820          m=is(imade)
 821          diff=dij(n,m)
 822          if(diff.lt.diffmn) diffmn=diff
 823         enddo
 824         if(diffmn.gt.difmax) difmax=diffmn
 825         adiff(iuse)=diffmn
 826         itag(iuse)=n
 827         avedif=avedif+diffmn
 828        endif
 829       enddo
 830
 831       avedif=avedif/iuse
 832 c     avedif=(avedif+difmax)/2
 833       emax=-9.d190
 834       do i=1,iuse
 835        if(adiff(i).ge.avedif) then
 836         itagi=itag(i)
 837         benei=bene(itagi)
 838         if(benei.gt.emax) then
 839          emax=benei
 840          is(imade1)=itagi
 841         endif
 842        endif
 843       enddo
 844
 845       if(ik.eq.0) iuse=iuse-1
 846
 847       return
 848       end
 849 c---------------------------------
 850       subroutine select_iseed_min(imade1,ik)
 851       implicit real*8 (a-h,o-z)
 852       include 'DIMENSIONS'
 853       include 'COMMON.CSA'
 854       include 'COMMON.BANK'
 855       dimension itag(mxio),adiff(mxio)
 856
 857       iuse=0
 858       avedif=0.d0
 859       difmax=0.d0
 860       do n=1,nbank
 861        if(ibank(n).eq.ik) then
 862         iuse=iuse+1
 863         diffmn=9.d190
 864         do imade=1,imade1-1
 865 c        m=nbank+imade
 866 c        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 867          m=is(imade)
 868          diff=dij(n,m)
 869          if(diff.lt.diffmn) diffmn=diff
 870         enddo
 871         if(diffmn.gt.difmax) difmax=diffmn
 872         adiff(iuse)=diffmn
 873         itag(iuse)=n
 874         avedif=avedif+diffmn
 875        endif
 876       enddo
 877
 878       avedif=avedif/iuse
 879 c     avedif=(avedif+difmax)/2
 880       emin=9.d190
 881       do i=1,iuse
 882 c      print *,"i, adiff(i),avedif : ",i,adiff(i),avedif
 883        if(adiff(i).ge.avedif) then
 884         itagi=itag(i)
 885         benei=bene(itagi)
 886 c       print *,"i, benei,emin : ",i,benei,emin
 887         if(benei.lt.emin) then
 888          emin=benei
 889          is(imade1)=itagi
 890         endif
 891        endif
 892       enddo
 893
 894       if(ik.eq.0) iuse=iuse-1
 895
 896 c     print *, "exiting select_iseed_min",is(imade1)
 897
 898       return
 899       end
 900 c---------------------------------
 901       subroutine select_iseed_far(imade1,ik)
 902       implicit real*8 (a-h,o-z)
 903       include 'DIMENSIONS'
 904       include 'COMMON.CSA'
 905       include 'COMMON.BANK'
 906
 907       dmax=-9.d190
 908       do n=1,nbank
 909        if(ibank(n).eq.ik) then
 910         diffmn=9.d190
 911         do imade=1,imade1-1
 912 c        m=nbank+imade
 913 c        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
 914          m=is(imade)
 915          diff=dij(n,m)
 916          if(diff.lt.diffmn) diffmn=diff
 917         enddo
 918        endif
 919        if(diffmn.gt.dmax) then
 920         dmax=diffmn
 921         is(imade1)=n
 922        endif
 923       enddo
 924
 925       return
 926       end
 927 c---------------------------------
 928       subroutine find_min
 929       implicit real*8 (a-h,o-z)
 930       include 'DIMENSIONS'
 931       include 'COMMON.CSA'
 932       include 'COMMON.BANK'
 933
 934       ebmin=9.d190
 935
 936       do i=1,nbank
 937        benei=bene(i)
 938        if(benei.lt.ebmin) then
 939         ebmin=benei
 940         ibmin=i
 941        endif
 942       enddo
 943
 944       return
 945       end
 946 c---------------------------------
 947       subroutine write_csa_pdb(var,ene,nft,ik,iw_pdb)
 948       implicit real*8 (a-h,o-z)
 949       include 'DIMENSIONS'
 950       include 'COMMON.CSA'
 951       include 'COMMON.BANK'
 952       include 'COMMON.VAR'
 953       include 'COMMON.IOUNITS'
 954       include 'COMMON.MINIM'
 955       include 'COMMON.SETUP'
 956       include 'COMMON.GEO'
 957       include 'COMMON.CHAIN'
 958       include 'COMMON.LOCAL'
 959       include 'COMMON.INTERACT'
 960       include 'COMMON.NAMES'
 961       include 'COMMON.SBRIDGE'
 962       integer lenpre,lenpot,ilen
 963       external ilen
 964       dimension var(maxvar)
 965       character*50 titelloc
 966       character*3 zahl
 967
 968       nmin_csa=nmin_csa+1
 969       if(ene.lt.eglob_csa) then
 970         eglob_csa=ene
 971         nglob_csa=nglob_csa+1
 972         call numstr(nglob_csa,zahl)
 973
 974         call var_to_geom(nvar,var)
 975         call chainbuild
 976         call secondary2(.false.)
 977
 978         lenpre=ilen(prefix)
 979         open(icsa_pdb,file=prefix(:lenpre)//'@'//zahl//'.pdb')
 980
 981         if (iw_pdb.eq.1) then
 982           write(titelloc,'(a2,i3,a3,i9,a3,i6)')
 983      &    'GM',nglob_csa,' e ',nft,' m ',nmin_csa
 984         else
 985           write(titelloc,'(a2,i3,a3,i9,a3,i6,a5,f5.2,a5,f5.1)')
 986      &   'GM',nglob_csa,' e ',nft,' m ',nmin_csa,' rms '
 987      &        ,rmsn(ik),' %NC ',pncn(ik)*100
 988         endif
 989         call pdbout(eglob_csa,titelloc,icsa_pdb)
 990         close(icsa_pdb)
 991       endif
 992
 993       return
 994       end
 995 c---------------------------------
 996       subroutine find_max
 997       implicit real*8 (a-h,o-z)
 998       include 'DIMENSIONS'
 999       include 'COMMON.CSA'
1000       include 'COMMON.BANK'
1001
1002       ebmax=-9.d190
1003
1004       do i=1,nbank
1005        benei=bene(i)
1006        if(benei.gt.ebmax) then
1007         ebmax=benei
1008         ibmax=i
1009        endif
1010       enddo
1011
1012       return
1013       end
1014 c---------------------------------
1015       subroutine get_diff
1016       implicit real*8 (a-h,o-z)
1017       include 'DIMENSIONS'
1018       include 'COMMON.CSA'
1019       include 'COMMON.BANK'
1020
1021       tdiff=0.d0
1022       difmin=9.d190
1023       do i1=1,nbank-1
1024        do i2=i1+1,nbank
1025         if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
1026          call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
1027          dij(i1,i2)=diff
1028          dij(i2,i1)=diff
1029         else
1030          diff=dij(i1,i2)
1031         endif
1032         tdiff=tdiff+diff
1033         if(diff.lt.difmin) difmin=diff
1034        enddo
1035        dij(i1,i1)=0.0
1036       enddo
1037
1038       do i=1,nbank
1039        jbank(i)=1
1040       enddo
1041
1042       avedif=tdiff/nbank/(nbank-1)*2
1043
1044       return
1045       end
1046 c---------------------------------
1047       subroutine estimate_cutdif(adif,xct,cutdifr)
1048       implicit real*8 (a-h,o-z)
1049       include 'DIMENSIONS'
1050       include 'COMMON.CSA'
1051       include 'COMMON.BANK'
1052
1053       ctdif1=adif/cut2
1054
1055       exponent = cutdifr*cut1/adif
1056       exponent = dlog(exponent)/dlog(xct)
1057
1058       nexp=exponent+0.25
1059       cutdif= adif/cut1*xct**nexp
1060       if(cutdif.lt.ctdif1) cutdif=ctdif1
1061
1062       return
1063       end
1064 c---------------------------------
1065       subroutine get_is_max(idum,ifar,n,imade,k)
1066       implicit real*8 (a-h,o-z)
1067       include 'DIMENSIONS'
1068       include 'COMMON.CSA'
1069       include 'COMMON.BANK'
1070       double precision emax
1071
1072       do i=imade+1,n
1073        emax=-9.d190
1074        do j=1,nbank
1075         if(ibank(j).eq.k .and. bene(j).gt.emax) then
1076            emax=bene(j)
1077            is(i)=j
1078         endif
1079        enddo
1080        call save_is(i)
1081       enddo
1082
1083       return
1084       end