source/unres/src_MD/src/bank.F

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