wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / 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