added source code
[unres.git] / 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