added v4.0 sources
[unres4.git] / source / unres / CSA.f90
1       module csa
2 !-----------------------------------------------------------------------------
3       use io_units
4       use names
5       use math
6       use random
7       use geometry_data, only: nres,rad2deg
8       use geometry
9       use energy
10       use MPI_
11       use csa_data
12       use compare
13       use io_config
14
15       implicit none
16 !-----------------------------------------------------------------------------
17 ! commom.bank
18 !      common/varin/
19 !      real(kind=8),dimension(:,:,:,:),allocatable :: dihang_in !(mxang,maxres,mxch,mxio)
20       integer,dimension(:),allocatable :: nss_in !(mxio)
21       integer,dimension(:,:),allocatable :: iss_in,jss_in !(maxss,mxio)
22 !      common/minvar/
23       real(kind=8),dimension(:,:,:,:),allocatable :: dihang !(mxang,maxres,mxch,mxio)
24       real(kind=8),dimension(:),allocatable :: etot!,rmsn,pncn !(mxio)
25       integer,dimension(:),allocatable :: nss_out !(mxio)
26       integer,dimension(:,:),allocatable ::iss_out,jss_out !(maxss,mxio)
27 !      common/bank/
28 !      real(kind=8),dimension(:,:,:,:),allocatable :: bvar !(mxang,maxres,mxch,mxio)
29 !      real(kind=8),dimension(:),allocatable :: bene,rene,&
30 !       brmsn,rrmsn,bpncn,rpncn !(mxio)
31       integer,dimension(:),allocatable :: is,jbank !(mxio)
32       real(kind=8) :: avedif,difmin,ebmin,ebmax,ebmaxt!,&
33 !       dele,difcut,cutdif,rmscut,pnccut
34       real(kind=8),dimension(:,:),allocatable :: dij !(mxio,mxio)
35 !      common/bank_disulfid/               
36 !      common/mvstat/
37       integer,dimension(:),allocatable :: movenx,movernx !(mxio)
38       integer,dimension(:,:),allocatable :: nstatnx,nstatnx_tot !(0:mxmv,3)
39       integer,dimension(:,:),allocatable :: indb !(mxio,9)
40       integer,dimension(:,:),allocatable :: parent !(3,mxio)
41 !      common/send2/
42       integer,dimension(:),allocatable :: isend2 !(mxio)
43       integer,dimension(:,:),allocatable :: iff_in !(maxres,mxio2)
44       integer,dimension(:,:,:,:),allocatable :: dihang_in2 !(mxang,maxres,mxch,mxio2)
45       integer,dimension(:,:),allocatable :: idata !(5,mxio)
46 !-----------------------------------------------------------------------------
47 ! common.csa
48 !      integer :: irestart,ndiff
49 !      common/alphaa/
50       integer,dimension(:),allocatable :: ngroup !(mxgr)
51       integer,dimension(:,:,:),allocatable :: igroup !(3,mxang,mxgr)
52       integer :: ntotgr!,numch
53 !      common/csa_input/
54 !      common/dih_control/
55 !      real(kind=8) :: rdih_bias
56       logical :: ldih_bias
57 !-----------------------------------------------------------------------------
58 ! common.distfit
59 !      COMMON /frag/
60       integer,dimension(:,:),allocatable :: bvar_frag !(mxio,6)
61       integer,dimension(:,:),allocatable :: hvar_frag,lvar_frag,svar_frag !(mxio,3)
62       integer,dimension(:,:),allocatable :: avar_frag !(mxio,5)
63 !-----------------------------------------------------------------------------
64 ! commom.hairpin
65 !      common /spinka/
66       integer :: nharp_tot
67       integer,dimension(:),allocatable :: nharp_seed,nharp_use  !(max_seed)
68       integer,dimension(:,:,:),allocatable :: iharp_seed        !(4,maxres/3,max_seed)
69       integer,dimension(:,:,:),allocatable :: iharp_use         !(0:4,maxres/3,max_seed)
70 !-----------------------------------------------------------------------------
71 ! Maximum number of moves (n1-n8)
72       integer,parameter :: mxmv=18
73 !-----------------------------------------------------------------------------
74 !
75 !
76 !-----------------------------------------------------------------------------
77       contains
78 !-----------------------------------------------------------------------------
79 ! bank.F
80 !-----------------------------------------------------------------------------
81       subroutine refresh_bank(ntrial)
82
83 !      implicit real*8 (a-h,o-z)
84 !      include 'DIMENSIONS'
85       include 'mpif.h'
86 !      include 'COMMON.CSA'
87 !      include 'COMMON.BANK'
88 !      include 'COMMON.IOUNITS'
89 !      include 'COMMON.CHAIN'
90 !      include 'COMMON.VAR'
91 !      include 'COMMON.CONTROL'
92       character :: chacc
93       integer :: iaccn,ntrial
94       real(kind=8) :: l_diff(mxio),denep
95       integer :: i,j,n,m,i1,idmin
96       real(kind=8) :: del_ene
97
98       do i=0,mxmv
99         do j=1,3
100           nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j)
101           nstatnx(i,j)=0
102         enddo
103       enddo
104
105 ! loop over all newly obtained conformations
106       do n=1,ntrial
107        chacc=' '
108        iaccn=0
109        nstatnx(movernx(n),1)=nstatnx(movernx(n),1)+1
110 !ccccccccccccccccccccccccccccccccccccccccccc
111 !jlee
112        if(iref.ne.0) then
113         if(rmsn(n).gt.rmscut.or.pncn(n).lt.pnccut) goto 100
114        endif
115 !jlee
116        if(etot(n).gt.ebmax) goto 100
117 ! Find the conformation closest to the conformation n in the bank
118        difmin=9.d9
119        do m=1,nbank
120         call get_diff12(dihang(1,1,1,n),bvar(1,1,1,m),l_diff(m))
121         if(l_diff(m).lt.difmin) then
122          difmin=l_diff(m)
123          idmin=m
124         endif
125        enddo
126
127        if(difmin.lt.cutdif) then
128 ! n is redundant to idmin
129         if(etot(n).lt.bene(idmin)) then
130          if(etot(n).lt.bene(idmin)-0.01d0) then
131           ibank(idmin)=0
132           jbank(idmin)=0
133          endif
134          denep=bene(idmin)-etot(n)
135          call replace_bvar(idmin,n)
136 !rc Update dij
137          do i1=1,nbank
138            if (i1.ne.idmin) then
139             dij(i1,idmin)=l_diff(i1)
140             dij(idmin,i1)=l_diff(i1)
141            endif
142          enddo
143          chacc='c'
144          iaccn=idmin
145          nstatnx(movernx(n),2)=nstatnx(movernx(n),2)+1
146          if(idmin.eq.ibmax) call find_max
147         endif
148        else
149 ! got new conformation
150         del_ene=0.0d0
151         if(ebmax-ebmin.gt.del_ene) then
152          denep=ebmax-etot(n)
153          call replace_bvar(ibmax,n)
154 !rc Update dij
155          do i1=1,nbank
156            if (i1.ne.ibmax) then
157             dij(i1,ibmax)=l_diff(i1)
158             dij(ibmax,i1)=l_diff(i1)
159            endif
160          enddo
161          chacc='f'
162          iaccn=ibmax
163          nstatnx(movernx(n),3)=nstatnx(movernx(n),3)+1
164          ibank(ibmax)=0
165          jbank(ibmax)=0
166          call find_max
167         else
168          if(del_ene.lt.0.0001) then
169           write (iout,*) 'ERROR in refresh_bank: '
170           write (iout,*) 'ebmax: ',ebmax
171           write (iout,*) 'ebmin: ',ebmin
172           write (iout,*) 'del_ene: ',del_ene
173 !rc          call mpi_abort(mpi_comm_world,ierror,ierrcode)
174          endif
175 !jp nbmax is never defined so condition below is always false
176 !         if(nbank.lt.nbmax) then
177 !          nbank=nbank+1
178 !          call replace_bvar(nbank,n)
179 !          ibank(nbank)=0
180 !          jbank(nbank)=0
181 !         else
182           call replace_bvar(ibmax,n)
183           ibank(ibmax)=0
184           jbank(ibmax)=0
185           call find_max
186 !         endif
187         endif
188        endif
189 !ccccccccccccccccccccccccccccccccccccccccccc
190   100 continue
191        if (iaccn.eq.0) then
192         if (iref.eq.0) then 
193          write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5)') &
194           indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
195           indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9)
196         else
197       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5,0pf4.1,a5,f3.0)') &
198           indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
199           indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
200           ' rms ',rmsn(n),' %NC ',pncn(n)*100
201         endif
202        else
203         if (iref.eq.0) then
204         write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,1x,a1,i4,0pf8.1,0pf8.1)') &
205           indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
206           indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
207           chacc,iaccn,difmin,denep
208         else
209          write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5,0pf4.1,&
210                                     a5,f3.0,1x,a1,i4,0pf8.1,0pf8.1)') &
211           indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',&
212           indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),&
213           ' rms ',rmsn(n),' %NC ',pncn(n)*100,&
214            chacc,iaccn,difmin,denep
215         endif
216        endif
217       enddo
218 ! end of loop over all newly obtained conformations
219       do i=0,mxmv
220         if(nstatnx(i,1).ne.0) then
221          if (i.le.9) then
222          write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
223                 '## N',i,' total=',nstatnx(i,1),&
224                 ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
225                 ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1)
226          else
227          write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
228                 '##N',i,' total=',nstatnx(i,1),&
229                 ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
230                 ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1)
231          endif
232         else
233          if (i.le.9) then        
234          write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
235                 '## N',i,' total=',nstatnx(i,1),&
236                 ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
237                 ' %acc',0.0
238          else
239          write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') &
240                 '##N',i,' total=',nstatnx(i,1),&
241                 ' close=',nstatnx(i,2),' far=',nstatnx(i,3),&
242                 ' %acc',0.0
243          endif
244         endif
245       enddo
246       call flush(iout)
247 !rc Update dij
248 !rc moved up, saves some get_diff12 calls 
249 !rc
250 !rc      do i1=1,nbank-1
251 !rc       do i2=i1+1,nbank
252 !rc        if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
253 !rc         call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
254 !rc         dij(i1,i2)=diff
255 !rc         dij(i2,i1)=diff
256 !rc        endif
257 !rc       enddo
258 !rc      enddo
259
260       do i=1,nbank
261        jbank(i)=1
262       enddo
263
264       return
265       end subroutine refresh_bank
266 !-----------------------------------------------------------------------------
267       subroutine replace_bvar(iold,inew)
268
269       use control_data, only: vdisulf
270       use energy_data, only: ns,iss
271 !      implicit real*8 (a-h,o-z)
272 !      include 'DIMENSIONS'
273       include 'mpif.h'
274 !      include 'COMMON.IOUNITS'
275 !      include 'COMMON.CSA'
276 !      include 'COMMON.BANK'
277 !      include 'COMMON.CHAIN'
278 !      include 'COMMON.CONTROL'
279 !      include 'COMMON.SBRIDGE'
280       integer :: iold,inew,ierror,ierrcode,i,j,k
281
282       if (iold.gt.mxio .or. iold.lt.1 .or. inew.gt.mxio .or. inew.lt.1) &
283         then
284         write (iout,*) 'Dimension ERROR in REPLACE_BVAR: IOLD',iold,&
285         ' INEW',inew
286         call mpi_abort(mpi_comm_world,ierror,ierrcode)
287       endif
288       do k=1,numch
289        do j=2,nres-1
290         do i=1,4
291          bvar(i,j,k,iold)=dihang(i,j,k,inew)
292         enddo
293        enddo
294       enddo
295       bene(iold)=etot(inew)
296       brmsn(iold)=rmsn(inew)
297       bpncn(iold)=pncn(inew)
298
299       if(bene(iold).lt.ebmin) then
300         ebmin=bene(iold)
301         ibmin=iold
302       endif
303
304       if(vdisulf) then
305         bvar_nss(iold)=nss_out(inew)
306 !d        write(iout,*) 'SS BANK',iold,bvar_nss(iold)
307         do i=1,bvar_nss(iold)
308           bvar_ss(1,i,iold)=iss_out(i,inew)
309           bvar_ss(2,i,iold)=jss_out(i,inew)
310 !d          write(iout,*) 'SS',bvar_ss(1,i,iold)-nres,
311 !d     &          bvar_ss(2,i,iold)-nres
312         enddo
313
314         bvar_ns(iold)=ns-2*bvar_nss(iold)
315 !d        write(iout,*) 'CYS #free ', bvar_ns(iold)
316            k=0
317            do i=1,ns
318              j=1
319              do while( iss(i).ne.iss_out(j,inew)-nres .and. &
320                        iss(i).ne.jss_out(j,inew)-nres .and. &
321                        j.le.nss_out(inew))
322                 j=j+1 
323              enddo
324              if (j.gt.nss_out(inew)) then            
325                k=k+1   
326                bvar_s(k,iold)=iss(i)
327              endif
328            enddo
329 !d         write(iout,*) 'CYS free',(bvar_s(k,iold),k=1,bvar_ns(iold))
330       endif
331
332       return
333       end subroutine replace_bvar
334 !-----------------------------------------------------------------------------
335       subroutine save_is(ind)
336
337 !      implicit real*8 (a-h,o-z)
338 !      include 'DIMENSIONS'
339       include 'mpif.h'
340 !      include 'COMMON.IOUNITS'
341 !      include 'COMMON.CSA'
342 !      include 'COMMON.BANK'
343 !      include 'COMMON.CHAIN'
344       integer :: ind,i,j,k,index,ierror,ierrcode
345
346       index=nbank+ind
347 !     print *, "nbank,ind,index,is(ind) ",nbank,ind,index,is(ind)
348       if (index.gt.mxio .or. index.lt.1 .or. &
349         is(ind).gt.mxio .or. is(ind).lt.1) then
350         write (iout,*) 'Dimension ERROR in SAVE_IS: INDEX',index,&
351         ' IND',ind,' IS',is(ind)
352         call mpi_abort(mpi_comm_world,ierror,ierrcode)
353       endif
354       do k=1,numch
355        do j=2,nres-1
356         do i=1,4
357          bvar(i,j,k,index)=bvar(i,j,k,is(ind))
358         enddo
359        enddo
360       enddo
361       bene(index)=bene(is(ind))
362       ibank(is(ind))=1
363
364       return
365       end subroutine save_is
366 !-----------------------------------------------------------------------------
367       subroutine select_is(n,ifar,idum)
368
369 !      implicit real*8 (a-h,o-z)
370 !      include 'DIMENSIONS'
371 !      include 'COMMON.CSA'
372 !      include 'COMMON.BANK'
373       integer,dimension(mxio) :: itag
374       real(kind=8),dimension(mxio) :: adiff
375       integer :: n,ifar,idum,i,iusesv,imade
376
377       iuse=0
378       do i=1,nbank
379        if(ibank(i).eq.0) then
380         iuse=iuse+1
381         itag(iuse)=i
382        endif
383       enddo
384       iusesv=iuse
385
386       if(iuse.eq.0) then
387        icycle=icycle+1
388        do i=1,nbank
389         if(ibank(i).eq.2) then
390          ibank(i)=1
391         else
392          ibank(i)=0
393         endif
394        enddo
395        imade=0
396        call get_is(idum,ifar,n,imade,0)
397 !test3       call get_is_max(idum,ifar,n,imade,0)
398       else if(iuse.eq.n) then
399        do i=1,iuse
400         is(i)=itag(i)
401         call save_is(i)
402        enddo
403       else if(iuse.lt.n) then
404 !      if(icycle.eq.0) then
405 !       do i=1,n
406 !        ind=mod(i-1,iuse)+1
407 !        is(i)=itag(ind)
408 !        call save_is(i)
409 !       enddo
410 !      else
411 !      endif
412        do i=1,iuse
413         is(i)=itag(i)
414         call save_is(i)
415        enddo
416        imade=iuse
417 !      call get_is_ran(idum,n,imade,1)
418        call get_is(idum,ifar,n,imade,1)
419 !test3       call get_is_max(idum,ifar,n,imade,1)
420 !      if(iusesv.le.n/10) then
421        if(iusesv.le.0) then
422         icycle=icycle+1
423         do i=1,nbank
424 !        if(ibank(i).eq.2) then
425 !         ibank(i)=1
426          if(ibank(i).ge.2) then
427           ibank(i)=ibank(i)-1
428          else
429           ibank(i)=0
430          endif
431         enddo
432        endif
433       else
434        imade=0
435        call get_is(idum,ifar,n,imade,0)
436 !test3       call get_is_max(idum,ifar,n,imade,0)
437       endif
438       iuse=iusesv
439
440       return
441       end subroutine select_is
442 !-----------------------------------------------------------------------------
443       subroutine get_is_ran(idum,n,imade,k)
444
445 !      implicit real*8 (a-h,o-z)
446 !      include 'DIMENSIONS'
447 !      include 'COMMON.CSA'
448 !      include 'COMMON.BANK'
449 !      real(kind=4) :: ran1,ran2
450       integer,dimension(mxio) :: itag
451       real(kind=8),dimension(mxio) :: adiff
452       integer :: idum,n,imade,k,j,i,iran
453
454       do j=imade+1,n
455        iuse=0
456        do i=1,nbank
457         if(ibank(i).eq.k) then
458          iuse=iuse+1
459          itag(iuse)=i
460         endif
461        enddo
462        iran=iuse*  ran1(idum)+1
463        is(j)=itag(iran)
464        call save_is(j)
465       enddo
466
467       return
468       end subroutine get_is_ran
469 !-----------------------------------------------------------------------------
470       subroutine get_is(idum,ifar,n,imade,k)
471
472 !      implicit real*8 (a-h,o-z)
473 !      include 'DIMENSIONS'
474 !      include 'COMMON.CSA'
475 !      include 'COMMON.BANK'
476 !      real(kind=4) :: ran1,ran2
477       integer,dimension(mxio) :: itag
478       real(kind=8),dimension(mxio) :: adiff
479       integer :: idum,ifar,n,imade,k,i,iran
480
481       iuse=0
482       do i=1,nbank
483        if(ibank(i).eq.k) then
484         iuse=iuse+1
485         itag(iuse)=i
486        endif
487       enddo
488       iran=iuse*  ran1(idum)+1
489       imade=imade+1
490       is(imade)=itag(iran)
491       call save_is(imade)
492
493       do i=imade+1,ifar-1
494        if(icycle.eq.-1) then
495         call select_iseed_max(i,k)
496        else
497         call select_iseed_min(i,k)
498 !test4  call select_iseed_max(i,k)
499        endif
500        call save_is(i)
501       enddo
502
503       do i=ifar,n
504        call select_iseed_far(i,k)
505        call save_is(i)
506       enddo
507
508       return
509       end subroutine get_is
510 !-----------------------------------------------------------------------------
511       subroutine select_iseed_max(imade1,ik)
512
513 !      implicit real*8 (a-h,o-z)
514 !      include 'DIMENSIONS'
515 !      include 'COMMON.CSA'
516 !      include 'COMMON.BANK'
517       integer,dimension(mxio) :: itag
518       real(kind=8),dimension(mxio) :: adiff
519       integer :: imade1,ik,i,n,imade,m,itagi
520       real(kind=8) :: difmax,diff,emax,benei,diffmn
521
522       iuse=0
523       avedif=0.d0
524       difmax=0.d0
525       do n=1,nbank
526        if(ibank(n).eq.ik) then
527         iuse=iuse+1
528         diffmn=9.d190
529         do imade=1,imade1-1
530 !        m=nbank+imade
531 !        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
532          m=is(imade)
533          diff=dij(n,m)
534          if(diff.lt.diffmn) diffmn=diff
535         enddo
536         if(diffmn.gt.difmax) difmax=diffmn
537         adiff(iuse)=diffmn
538         itag(iuse)=n
539         avedif=avedif+diffmn
540        endif
541       enddo
542
543       avedif=avedif/iuse
544 !     avedif=(avedif+difmax)/2
545       emax=-9.d190
546       do i=1,iuse
547        if(adiff(i).ge.avedif) then
548         itagi=itag(i)
549         benei=bene(itagi)
550         if(benei.gt.emax) then
551          emax=benei
552          is(imade1)=itagi  
553         endif
554        endif
555       enddo
556
557       if(ik.eq.0) iuse=iuse-1
558
559       return
560       end subroutine select_iseed_max
561 !-----------------------------------------------------------------------------
562       subroutine select_iseed_min(imade1,ik)
563
564 !      implicit real*8 (a-h,o-z)
565 !      include 'DIMENSIONS'
566 !      include 'COMMON.CSA'
567 !      include 'COMMON.BANK'
568       integer,dimension(mxio) :: itag
569       real(kind=8),dimension(mxio) :: adiff
570       integer :: imade1,ik,n,imade,m,i,itagi
571       real(kind=8) :: difmax,diff,diffmn,emin,benei
572
573       iuse=0
574       avedif=0.d0
575       difmax=0.d0
576       do n=1,nbank
577        if(ibank(n).eq.ik) then
578         iuse=iuse+1
579         diffmn=9.d190
580         do imade=1,imade1-1
581 !        m=nbank+imade
582 !        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
583          m=is(imade)
584          diff=dij(n,m)
585          if(diff.lt.diffmn) diffmn=diff
586         enddo
587         if(diffmn.gt.difmax) difmax=diffmn
588         adiff(iuse)=diffmn
589         itag(iuse)=n
590         avedif=avedif+diffmn
591        endif
592       enddo
593
594       avedif=avedif/iuse
595 !     avedif=(avedif+difmax)/2
596       emin=9.d190
597       do i=1,iuse
598 !      print *,"i, adiff(i),avedif : ",i,adiff(i),avedif
599        if(adiff(i).ge.avedif) then
600         itagi=itag(i)
601         benei=bene(itagi)
602 !       print *,"i, benei,emin : ",i,benei,emin
603         if(benei.lt.emin) then
604          emin=benei
605          is(imade1)=itagi  
606         endif
607        endif
608       enddo
609
610       if(ik.eq.0) iuse=iuse-1
611
612 !     print *, "exiting select_iseed_min",is(imade1)
613
614       return
615       end subroutine select_iseed_min
616 !-----------------------------------------------------------------------------
617       subroutine select_iseed_far(imade1,ik)
618
619 !      implicit real*8 (a-h,o-z)
620 !      include 'DIMENSIONS'
621 !      include 'COMMON.CSA'
622 !      include 'COMMON.BANK'
623       integer :: imade1,ik,n,imade,m
624       real(kind=8) :: dmax,diffmn,diff
625
626       dmax=-9.d190
627       do n=1,nbank
628        if(ibank(n).eq.ik) then
629         diffmn=9.d190
630         do imade=1,imade1-1
631 !        m=nbank+imade
632 !        call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff)
633          m=is(imade)
634          diff=dij(n,m)
635          if(diff.lt.diffmn) diffmn=diff
636         enddo
637        endif
638        if(diffmn.gt.dmax) then
639         dmax=diffmn
640         is(imade1)=n  
641        endif
642       enddo
643
644       return
645       end subroutine select_iseed_far
646 !-----------------------------------------------------------------------------
647       subroutine find_min
648
649 !      implicit real*8 (a-h,o-z)
650 !      include 'DIMENSIONS'
651 !      include 'COMMON.CSA'
652 !      include 'COMMON.BANK'
653       integer :: i
654       real(kind=8) :: benei
655
656       ebmin=9.d190
657
658       do i=1,nbank
659        benei=bene(i)
660        if(benei.lt.ebmin) then
661         ebmin=benei
662         ibmin=i
663        endif   
664       enddo    
665
666       return
667       end subroutine find_min
668 !-----------------------------------------------------------------------------
669       subroutine find_max
670
671 !      implicit real*8 (a-h,o-z)
672 !      include 'DIMENSIONS'
673 !      include 'COMMON.CSA'
674 !      include 'COMMON.BANK'
675       integer :: i
676       real(kind=8) :: benei
677
678       ebmax=-9.d190
679
680       do i=1,nbank
681        benei=bene(i)
682        if(benei.gt.ebmax) then
683         ebmax=benei
684         ibmax=i
685        endif
686       enddo
687
688       return
689       end subroutine find_max
690 !-----------------------------------------------------------------------------
691       subroutine get_diff
692
693 !      implicit real*8 (a-h,o-z)
694 !      include 'DIMENSIONS'
695 !      include 'COMMON.CSA'
696 !      include 'COMMON.BANK'
697       integer :: i,i1,i2
698       real(kind=8) :: tdiff,difmin,diff
699
700       tdiff=0.d0
701       difmin=9.d190
702       do i1=1,nbank-1
703        do i2=i1+1,nbank
704         if(jbank(i1).eq.0.or.jbank(i2).eq.0) then
705          call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff)
706          dij(i1,i2)=diff
707          dij(i2,i1)=diff
708         else
709          diff=dij(i1,i2)
710         endif
711         tdiff=tdiff+diff
712         if(diff.lt.difmin) difmin=diff
713        enddo
714        dij(i1,i1)=0.0
715       enddo
716
717       do i=1,nbank
718        jbank(i)=1
719       enddo
720
721       avedif=tdiff/nbank/(nbank-1)*2
722
723       return
724       end subroutine get_diff
725 !-----------------------------------------------------------------------------
726       subroutine estimate_cutdif(adif,xct,cutdifr)
727
728 !      implicit real*8 (a-h,o-z)
729 !      include 'DIMENSIONS'
730 !      include 'COMMON.CSA'
731 !      include 'COMMON.BANK'
732       integer :: nexp
733       real(kind=8) :: adif,xct,cutdifr,ctdif1,exponent
734
735       ctdif1=adif/cut2
736
737       exponent = cutdifr*cut1/adif
738       exponent = dlog(exponent)/dlog(xct)
739
740       nexp=exponent+0.25
741       cutdif= adif/cut1*xct**nexp
742       if(cutdif.lt.ctdif1) cutdif=ctdif1
743
744       return
745       end subroutine estimate_cutdif
746 !-----------------------------------------------------------------------------
747       subroutine get_is_max(idum,ifar,n,imade,k)
748
749 !      implicit real*8 (a-h,o-z)
750 !      include 'DIMENSIONS'
751 !      include 'COMMON.CSA'
752 !      include 'COMMON.BANK'
753       integer :: idum,ifar,n,imade,k,i,j
754       real(kind=8) :: emax
755
756       do i=imade+1,n
757        emax=-9.d190
758        do j=1,nbank
759         if(ibank(j).eq.k .and. bene(j).gt.emax) then
760            emax=bene(j)
761            is(i)=j
762         endif
763        enddo
764        call save_is(i)
765       enddo
766
767       return
768       end subroutine get_is_max
769 !-----------------------------------------------------------------------------
770 ! csa.f
771 !-----------------------------------------------------------------------------
772       subroutine make_array
773
774       use energy_data, only: itype
775 !      implicit real*8 (a-h,o-z)
776 !      include 'DIMENSIONS'
777 !      include 'COMMON.IOUNITS'
778 !      include 'COMMON.CHAIN'
779 !      include 'COMMON.INTERACT'
780 !      include 'COMMON.CSA'
781        integer :: k,j,i,indg
782 !cccccccccccccccccccccccc
783 ! Level-2: group
784 !cccccccccccccccccccccccc
785
786       indg=0
787       do k=1,numch
788 !cccccccccccccccccccccccccccccccccccccccc
789 ! Groups the THETAs and the GAMMAs
790        do j=2,nres-1
791          indg=indg+1
792          if (j.lt.nres-1) then
793            ngroup(indg)=2
794          else
795            ngroup(indg)=1
796          endif
797          do i=1,ngroup(indg)
798           igroup(1,i,indg)=i
799           igroup(2,i,indg)=j
800           igroup(3,i,indg)=k
801          enddo
802        enddo
803 !cccccccccccccccccccccccccccccccccccccccc
804       enddo
805 ! Groups the ALPHAs and the BETAs
806       do k=1,numch
807        do j=2,nres-1
808         if(itype(j).ne.10) then
809          indg=indg+1
810          ngroup(indg)=2
811          do i=1,ngroup(indg)
812           igroup(1,i,indg)=i+2
813           igroup(2,i,indg)=j
814           igroup(3,i,indg)=k
815          enddo
816         endif
817        enddo
818       enddo
819
820       ntotgr=indg
821       write(iout,*) 
822       write(iout,*) "# of groups: ",ntotgr
823       do i=1,ntotgr
824        write(iout,41) i,ngroup(i),((igroup(k,j,i),k=1,3),j=1,ngroup(i))
825       enddo
826 !      close(iout)
827
828    40 format(i3,3x,3i3)
829    41 format(2i3,3x,6(3i3,2x))
830
831       return
832       end subroutine make_array
833 !-----------------------------------------------------------------------------
834       subroutine make_ranvar(n,m,idum)
835
836       use geometry_data
837 !      implicit real*8 (a-h,o-z)
838 !      include 'DIMENSIONS'
839 !      include 'COMMON.IOUNITS'
840 !      include 'COMMON.CHAIN'
841 !      include 'COMMON.VAR'
842 !      include 'COMMON.BANK'
843       integer :: n,m,j,idum,itrial,jeden
844
845 ! al      m=0
846       print *,'HOHOHOHO Make_RanVar!!!!!',n,m
847       itrial=0
848       do while(m.lt.n .and. itrial.le.10000)
849         itrial=itrial+1
850         jeden=1
851         call gen_rand_conf(jeden,*10)
852 !        call intout
853         m=m+1
854         do j=2,nres-1
855           dihang_in(1,j,1,m)=theta(j+1)
856           dihang_in(2,j,1,m)=phi(j+2)
857           dihang_in(3,j,1,m)=alph(j)
858           dihang_in(4,j,1,m)=omeg(j)
859         enddo  
860         dihang_in(2,nres-1,1,m)=0.0d0
861         goto 20
862   10    write (iout,*) 'Failed to generate conformation #',m+1,&
863          ' itrial=',itrial
864   20    continue
865       enddo
866       print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial
867
868       return
869       end subroutine make_ranvar
870 !-----------------------------------------------------------------------------
871       subroutine make_ranvar_reg(n,idum)
872
873       use geometry_data
874 !      implicit real*8 (a-h,o-z)
875 !      include 'DIMENSIONS'
876 !      include 'COMMON.IOUNITS'
877 !      include 'COMMON.CHAIN'
878 !      include 'COMMON.VAR'
879 !      include 'COMMON.BANK'
880 !      include 'COMMON.GEO'
881       integer :: n,idum,j,m,itrial,jeden
882
883       m=0
884       print *,'HOHOHOHO Make_RanVar!!!!!'
885       itrial=0
886       do while(m.lt.n .and. itrial.le.10000)
887         itrial=itrial+1
888         jeden=1
889         call gen_rand_conf(jeden,*10)
890 !        call intout
891         m=m+1
892         do j=2,nres-1
893           dihang_in(1,j,1,m)=theta(j+1)
894           dihang_in(2,j,1,m)=phi(j+2)
895           dihang_in(3,j,1,m)=alph(j)
896           dihang_in(4,j,1,m)=omeg(j)
897           if(m.le.n*0.1) then
898            dihang_in(1,j,1,m)=90.0*deg2rad
899            dihang_in(2,j,1,m)=50.0*deg2rad
900           endif
901         enddo  
902         dihang_in(2,nres-1,1,m)=0.0d0
903         goto 20
904   10    write (iout,*) 'Failed to generate conformation #',m+1,&
905          ' itrial=',itrial
906   20    continue
907       enddo
908       print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial
909
910       return
911       end subroutine make_ranvar_reg
912 !-----------------------------------------------------------------------------
913 ! diff12.f
914 !-----------------------------------------------------------------------------
915       subroutine get_diff12(aarray,barray,diff)
916
917 !      implicit real*8 (a-h,o-z)
918 !      include 'DIMENSIONS'
919 !      include 'COMMON.CSA'
920 !      include 'COMMON.BANK'
921 !      include 'COMMON.CHAIN'
922 !      include 'COMMON.GEO'
923       integer :: k,j,i
924       real(kind=8),dimension(mxang,nres,mxch) :: aarray,barray  !(mxang,maxres,mxch)
925       real(kind=8) :: diff,dif
926
927       diff=0.d0
928       do k=1,numch
929        do j=2,nres-1
930 !       do i=1,4
931 !        do i=1,2
932         do i=1,ndiff
933          dif=rad2deg*dabs(aarray(i,j,k)-barray(i,j,k))
934          if(dif.gt.180.) dif=360.-dif
935          if (dif.gt.diffcut) diff=diff+dif
936         enddo
937        enddo
938       enddo
939
940       return
941       end subroutine get_diff12
942 !-----------------------------------------------------------------------------
943 ! indexx.f
944 !-----------------------------------------------------------------------------
945       subroutine indexx(n,arr,indx)
946
947 !      implicit real*8 (a-h,o-z)
948       INTEGER :: n,indx(n)
949       REAL(kind=8) :: arr(n)
950 !     PARAMETER (M=7,NSTACK=50)
951       integer,PARAMETER :: M=7,NSTACK=500
952       INTEGER :: i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
953       REAL(kind=8) :: a
954
955       do 11 j=1,n
956         indx(j)=j
957 11    continue
958       jstack=0
959       l=1
960       ir=n
961 1     if(ir-l.lt.M)then
962         do 13 j=l+1,ir
963           indxt=indx(j)
964           a=arr(indxt)
965           do 12 i=j-1,1,-1
966             if(arr(indx(i)).le.a)goto 2
967             indx(i+1)=indx(i)
968 12        continue
969           i=0
970 2         indx(i+1)=indxt
971 13      continue
972         if(jstack.eq.0)return
973         ir=istack(jstack)
974         l=istack(jstack-1)
975         jstack=jstack-2
976       else
977         k=(l+ir)/2
978         itemp=indx(k)
979         indx(k)=indx(l+1)
980         indx(l+1)=itemp
981         if(arr(indx(l+1)).gt.arr(indx(ir)))then
982           itemp=indx(l+1)
983           indx(l+1)=indx(ir)
984           indx(ir)=itemp
985         endif
986         if(arr(indx(l)).gt.arr(indx(ir)))then
987           itemp=indx(l)
988           indx(l)=indx(ir)
989           indx(ir)=itemp
990         endif
991         if(arr(indx(l+1)).gt.arr(indx(l)))then
992           itemp=indx(l+1)
993           indx(l+1)=indx(l)
994           indx(l)=itemp
995         endif
996         i=l+1
997         j=ir
998         indxt=indx(l)
999         a=arr(indxt)
1000 3       continue
1001           i=i+1
1002         if(arr(indx(i)).lt.a)goto 3
1003 4       continue
1004           j=j-1
1005         if(arr(indx(j)).gt.a)goto 4
1006         if(j.lt.i)goto 5
1007         itemp=indx(i)
1008         indx(i)=indx(j)
1009         indx(j)=itemp
1010         goto 3
1011 5       indx(l)=indx(j)
1012         indx(j)=indxt
1013         jstack=jstack+2
1014         if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx'
1015         if(ir-i+1.ge.j-l)then
1016           istack(jstack)=ir
1017           istack(jstack-1)=i
1018           ir=j-1
1019         else
1020           istack(jstack)=j-1
1021           istack(jstack-1)=l
1022           l=i
1023         endif
1024       endif
1025       goto 1
1026       end subroutine indexx
1027 !  (C) Copr. 1986-92 Numerical Recipes Software *11915aZ%.
1028 !-----------------------------------------------------------------------------
1029 ! minim_jlee.F
1030 !-----------------------------------------------------------------------------
1031       subroutine minim_jlee
1032
1033       use minim_data
1034       use MPI_data
1035       use energy_data
1036       use compare_data
1037       use control_data
1038       use geometry_data, only: nvar,nphi
1039       use geometry, only:dist
1040       use energy, only:fdum
1041       use control, only:init_int_table
1042       use minimm, only:sumsl,deflt
1043 !  controls minimization and sorting routines
1044 !      implicit real*8 (a-h,o-z)
1045 !      include 'DIMENSIONS'
1046 !      include 'COMMON.VAR'
1047 !      include 'COMMON.IOUNITS'
1048 !      include 'COMMON.MINIM'
1049 !      include 'COMMON.CONTROL'
1050       include 'mpif.h'
1051       integer,parameter :: liv=60
1052       integer :: lv
1053 !      external func,gradient!,fdum     !use minim & energy
1054 !      real(kind=4) :: ran1,ran2,ran3
1055 !      include 'COMMON.SETUP'
1056 !      include 'COMMON.GEO'
1057 !      include 'COMMON.FFIELD'
1058 !      include 'COMMON.SBRIDGE'
1059 !      include 'COMMON.DISTFIT'
1060 !      include 'COMMON.CHAIN'
1061       integer,dimension(mpi_status_size) :: muster
1062       real(kind=8),dimension(6*nres) :: var     !(maxvar) (maxvar=6*maxres)
1063       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: erg
1064       real(kind=8),dimension(6*nres) :: var2    !(maxvar) (maxvar=6*maxres)
1065       integer,dimension(nres) :: iffr   !(maxres)
1066       integer,dimension((nres-1)*(nres-2)/2) :: ihpbt,jhpbt     !(maxdim) (maxdim=(maxres-1)*(maxres-2)/2)
1067       real(kind=8),dimension(6*nres) :: d,garbage       !(maxvar) (maxvar=6*maxres)
1068 !el      real(kind=8),dimension(1:lv+1) :: v                 
1069       real(kind=8) :: energia(0:n_ene),time0s,time1s
1070       integer,dimension(9) :: indx
1071       integer,dimension(12) :: info
1072       integer,dimension(liv) :: iv                                               
1073       integer :: idum(1)
1074       real(kind=8) :: rdum(1)
1075       integer,dimension(2,12*nres) :: icont_    !(2,maxcont)(maxcont=12*maxres)
1076       logical :: fail   !check_var,
1077       integer :: iloop(2)
1078 !el      common /przechowalnia/ v
1079       integer :: i,j,ierr,n,nfun,nft_sc,nf,ierror,ierrcode
1080       real(kind=8) :: rad,eee,etot      !,fdum
1081 !el  from subroutine parmread
1082 ! Define the constants of the disulfide bridge
1083 ! Old arbitrary potential
1084       real(kind=8),parameter :: dbr=4.20D0
1085       real(kind=8),parameter :: fbr=3.30D0
1086 !-----------------
1087       lv=77+(6*nres)*(6*nres+17)/2      !77+maxvar*(maxvar+17)/2 (maxvar=6*maxres)
1088       data rad /1.745329252d-2/
1089 !  receive # of start
1090 !      print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun,
1091 !     &   ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf
1092       if (.not. allocated(v)) allocate(v(1:lv))
1093       nhpb0=nhpb
1094    10 continue
1095       time0s=MPI_WTIME()
1096 !      print *, 'MINIM_JLEE: ',me,' is waiting'
1097       call mpi_recv(info,12,mpi_integer,king,idint,CG_COMM,&
1098                     muster,ierr)
1099       time1s=MPI_WTIME()
1100       write (iout,'(a12,f10.4,a4)')'Waiting for ',time1s-time0s,' sec'
1101       call flush(iout)
1102        n=info(1)
1103 !      print *, 'MINIM_JLEE: ',me,' received: ',n
1104       
1105 !rc      if (ierr.ne.0) go to 100
1106 !  if # = 0, return
1107       if (n.eq.0) then 
1108         write (iout,*) 'Finishing minim_jlee - signal',n,' from master'
1109         call flush(iout)
1110         return
1111       endif
1112
1113       nfun=0
1114       IF (n.lt.0) THEN
1115        call mpi_recv(var,nvar,mpi_double_precision,&
1116                     king,idreal,CG_COMM,muster,ierr)
1117        call mpi_recv(iffr,nres,mpi_integer,&
1118                     king,idint,CG_COMM,muster,ierr)
1119        call mpi_recv(var2,nvar,mpi_double_precision,&
1120                     king,idreal,CG_COMM,muster,ierr)
1121       ELSE
1122 !  receive initial values of variables
1123        call mpi_recv(var,nvar,mpi_double_precision,&
1124                     king,idreal,CG_COMM,muster,ierr)
1125 !rc       if (ierr.ne.0) go to 100
1126       ENDIF
1127
1128       if(vdisulf.and.info(2).ne.-1) then
1129         if(info(4).ne.0)then
1130            call mpi_recv(ihpbt,info(4),mpi_integer,&
1131                     king,idint,CG_COMM,muster,ierr)
1132            call mpi_recv(jhpbt,info(4),mpi_integer,&
1133                     king,idint,CG_COMM,muster,ierr)
1134         endif
1135       endif  
1136
1137       IF (n.lt.0) THEN
1138         n=-n     
1139         nhpb=nhpb0
1140         link_start=1
1141         link_end=nhpb
1142         call init_int_table
1143         call contact_cp(var,var2,iffr,nfun,n)
1144       ENDIF
1145
1146       if(vdisulf.and.info(2).ne.-1) then
1147         nss=0
1148         if(info(4).ne.0)then
1149 !d          write(iout,*) 'SS=',info(4),'N=',info(1),'IT=',info(2)
1150           call var_to_geom(nvar,var)
1151           call chainbuild
1152           do i=1,info(4)
1153            if (dist(ihpbt(i),jhpbt(i)).lt.7.0) then
1154             nss=nss+1
1155             ihpb(nss)=ihpbt(i)
1156             jhpb(nss)=jhpbt(i)
1157 !d            write(iout,*) 'SS  mv=',info(3),
1158 !d     &         ihpb(nss)-nres,jhpb(nss)-nres,
1159 !d     &         dist(ihpb(nss),jhpb(nss))
1160             dhpb(nss)=dbr
1161             forcon(nss)=fbr
1162            else
1163 !d            write(iout,*) 'rm SS  mv=',info(3),
1164 !d     &         ihpbt(i)-nres,jhpbt(i)-nres,dist(ihpbt(i),jhpbt(i))
1165            endif
1166           enddo
1167         endif
1168         nhpb=nss
1169         link_start=1
1170         link_end=nhpb
1171         call init_int_table
1172       endif
1173
1174       if (info(3).eq.14) then
1175        write(iout,*) 'calling local_move',info(7),info(8)
1176        call local_move_init(.false.)
1177        call var_to_geom(nvar,var)
1178        call local_move(info(7),info(8),20d0,50d0)
1179        call geom_to_var(nvar,var)
1180       endif
1181
1182
1183       if (info(3).eq.16) then
1184        write(iout,*) 'calling beta_slide',info(7),info(8),&
1185                   info(10), info(11), info(12)
1186        call var_to_geom(nvar,var)
1187        call beta_slide(info(7),info(8),info(10),info(11),info(12), &
1188                                                            nfun,n)
1189        call geom_to_var(nvar,var)
1190       endif
1191
1192
1193       if (info(3).eq.17) then
1194        write(iout,*) 'calling beta_zip',info(7),info(8)
1195        call var_to_geom(nvar,var)
1196        call beta_zip(info(7),info(8),nfun,n)
1197        call geom_to_var(nvar,var)
1198       endif
1199
1200
1201 !rc overlap test
1202
1203       if (overlapsc) then 
1204
1205           call var_to_geom(nvar,var)
1206           call chainbuild
1207           call etotal(energia)
1208           nfun=nfun+1
1209           if (energia(1).eq.1.0d20) then  
1210             info(3)=-info(3) 
1211             write (iout,'(a,1pe14.5)')'#OVERLAP evdw=1d20',energia(1)
1212             call overlap_sc(fail)
1213             if(.not.fail) then
1214              call geom_to_var(nvar,var)
1215              call etotal(energia)
1216              nfun=nfun+1
1217              write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1)
1218             else
1219              v(10)=1.0d20
1220              iv(1)=-1
1221              goto 201
1222             endif
1223           endif
1224       endif 
1225
1226       if (searchsc) then 
1227           call var_to_geom(nvar,var)
1228           call sc_move(2,nres-1,1,10d0,nft_sc,etot)
1229           call geom_to_var(nvar,var)
1230 !d          write(iout,*) 'sc_move',nft_sc,etot
1231       endif 
1232
1233       if (check_var(var,info)) then 
1234            v(10)=1.0d21
1235            iv(1)=6
1236            goto 201
1237       endif
1238
1239
1240 !rc        
1241
1242 !      write (iout,*) 'MINIM_JLEE: Processor',me,' nvar',nvar
1243 !      write (iout,'(8f10.4)') (var(i),i=1,nvar)
1244 !      write (*,*) 'MINIM_JLEE: Processor',me,' received nvar',nvar
1245 !      write (*,'(8f10.4)') (var(i),i=1,nvar)
1246
1247        do i=1,nvar
1248          garbage(i)=var(i)
1249        enddo
1250
1251       call deflt(2,iv,liv,lv,v)                                         
1252 ! 12 means fresh start, dont call deflt                                 
1253       iv(1)=12                                                          
1254 ! max num of fun calls                                                  
1255       if (maxfun.eq.0) maxfun=500
1256       iv(17)=maxfun
1257 ! max num of iterations                                                 
1258       if (maxmin.eq.0) maxmin=1000
1259       iv(18)=maxmin
1260 ! controls output                                                       
1261       iv(19)=2                                                          
1262 ! selects output unit                                                   
1263 !d      iv(21)=iout                                                       
1264       iv(21)=0
1265 ! 1 means to print out result                                           
1266       iv(22)=0                                                          
1267 !d        iv(22)=1
1268 ! 1 means to print out summary stats                                    
1269       iv(23)=0                                                          
1270 ! 1 means to print initial x and d                                      
1271       iv(24)=0                                                          
1272
1273 !      if(me.eq.3.and.n.eq.255) then
1274 !       print *,' CHUJ: stoi'
1275 !       iv(21)=6
1276 !       iv(22)=1
1277 !       iv(23)=1
1278 !       iv(24)=1                                                          
1279 !      endif
1280
1281 ! min val for v(radfac) default is 0.1                                  
1282       v(24)=0.1D0                                                       
1283 ! max val for v(radfac) default is 4.0                                  
1284       v(25)=2.0D0                                                       
1285 !     v(25)=4.0D0                                                       
1286 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
1287 ! the sumsl default is 0.1                                              
1288       v(26)=0.1D0
1289 ! false conv if (act fnctn decrease) .lt. v(34)                         
1290 ! the sumsl default is 100*machep                                       
1291       v(34)=v(34)/100.0D0                                               
1292 ! absolute convergence                                                  
1293       if (tolf.eq.0.0D0) tolf=1.0D-4
1294       v(31)=tolf
1295 ! relative convergence                                                  
1296       if (rtolf.eq.0.0D0) rtolf=1.0D-4
1297       v(32)=rtolf
1298 ! controls initial step size                                            
1299        v(35)=1.0D-1                                                    
1300 ! large vals of d correspond to small components of step                
1301       do i=1,nphi
1302         d(i)=1.0D-1
1303       enddo
1304       do i=nphi+1,nvar
1305         d(i)=1.0D-1
1306       enddo
1307 !  minimize energy
1308 !      write (iout,*) 'Processor',me,' nvar',nvar
1309 !      write (iout,*) 'Variables BEFORE minimization:'
1310 !      write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
1311
1312 !      print *, 'MINIM_JLEE: ',me,' before SUMSL '
1313
1314       call func(nvar,var,nf,eee,idum,rdum,fdum)
1315       nfun=nfun+1
1316       if(eee.ge.1.0d20) then
1317 !       print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
1318 !       print *,' energy before SUMSL =',eee
1319 !       print *,' aborting local minimization'
1320        iv(1)=-1
1321        v(10)=eee
1322        go to 201
1323       endif
1324
1325 !t      time0s=MPI_WTIME()
1326       call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
1327 !t      write(iout,*) 'sumsl time=',MPI_WTIME()-time0s,iv(7),v(10)
1328 !      print *, 'MINIM_JLEE: ',me,' after SUMSL '
1329
1330 !  find which conformation was returned from sumsl
1331         nfun=nfun+iv(7)
1332 !      print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf,
1333 !     & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32)
1334 !        if (iv(1).ne.4 .or. nf.le.1) then
1335 !          write (*,*) 'Processor',me,' something bad in SUMSL',iv(1),nf
1336 !         write (*,*) 'Initial Variables'
1337 !         write (*,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
1338 !         write (*,*) 'Variables'
1339 !         write (*,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
1340 !         write (*,*) 'Vector d'
1341 !         write (*,'(8f10.4)') (d(i),i=1,nvar)
1342 !         write (iout,*) 'Processor',me,' something bad in SUMSL',
1343 !    &       iv(1),nf
1344 !         write (iout,*) 'Initial Variables'
1345 !         write (iout,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
1346 !         write (iout,*) 'Variables'
1347 !         write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
1348 !         write (iout,*) 'Vector d'
1349 !         write (iout,'(8f10.4)') (d(i),i=1,nvar)
1350 !        endif
1351 !        if (nf.lt.iv(6)-1) then
1352 !  recalculate intra- and interchain energies
1353 !         call func(nvar,var,nf,v(10),iv,v,fdum)
1354 !        else if (nf.eq.iv(6)-1) then
1355 !  regenerate conformation
1356 !         call var_to_geom(nvar,var)
1357 !         call chainbuild
1358 !        endif
1359 !  change origin and axes to standard ECEPP format
1360 !        call var_to_geom(nvar,var)
1361 !      write (iout,*) 'MINIM_JLEE after minim: Processor',me,' nvar',nvar
1362 !      write (iout,'(8f10.4)') (var(i),i=1,nvar)
1363 !      write (iout,*) 'Energy:',v(10)
1364 !  send back output
1365 !       print *, 'MINIM_JLEE: ',me,' minimized: ',n
1366   201  continue
1367         indx(1)=n
1368 ! return code: 6-gradient 9-number of ftn evaluation, etc
1369         indx(2)=iv(1)
1370 ! total # of ftn evaluations (for iwf=0, it includes all minimizations).
1371         indx(3)=nfun
1372         indx(4)=info(2)
1373         indx(5)=info(3)
1374         indx(6)=nss
1375         indx(7)=info(5)
1376         indx(8)=info(6)
1377         indx(9)=info(9)
1378         call mpi_send(indx,9,mpi_integer,king,idint,CG_COMM,&
1379                        ierr)
1380 !  send back energies
1381 ! al & cc
1382 ! calculate contact order
1383 #ifdef CO_BIAS
1384         call contact(.false.,ncont,icont_,co)
1385         erg(1)=v(10)-1.0d2*co
1386 #else
1387         erg(1)=v(10)
1388 #endif
1389         j=1
1390         call mpi_send(erg,j,mpi_double_precision,king,idreal,&
1391                        CG_COMM,ierr)
1392 #ifdef CO_BIAS
1393         call mpi_send(co,j,mpi_double_precision,king,idreal,&
1394                        CG_COMM,ierr)
1395 #endif
1396 !  send back values of variables
1397         call mpi_send(var,nvar,mpi_double_precision,&
1398                      king,idreal,CG_COMM,ierr)
1399 !        print * , 'MINIM_JLEE: Processor',me,' send erg and var '
1400
1401         if(vdisulf.and.info(2).ne.-1.and.nss.ne.0) then
1402 !d          call intout
1403 !d          call chainbuild
1404 !d          call etotal(energia(0))
1405 !d          etot=energia(0)
1406 !d          call enerprint(energia(0))
1407          call mpi_send(ihpb,nss,mpi_integer,&
1408                      king,idint,CG_COMM,ierr)        
1409          call mpi_send(jhpb,nss,mpi_integer,&
1410                      king,idint,CG_COMM,ierr)        
1411         endif
1412
1413         go to 10
1414   100 print *, ' error in receiving message from emperor', me
1415       call mpi_abort(mpi_comm_world,ierror,ierrcode)
1416       return
1417   200 print *, ' error in sending message to emperor'
1418       call mpi_abort(mpi_comm_world,ierror,ierrcode)
1419       return
1420   300 print *, ' error in communicating with emperor'
1421       call mpi_abort(mpi_comm_world,ierror,ierrcode)
1422       return
1423   956 format (' initial energy could not be calculated',41x)
1424   957 format (80x)
1425   965 format (' convergence code ',i2,' # of function calls ',&
1426         i4,' # of gradient calls ',i4,10x)
1427   975 format (' energy ',1p,e12.4,' scaled gradient ',e11.3,32x)
1428       end subroutine minim_jlee
1429 !-----------------------------------------------------------------------------
1430 ! newconf.f
1431 !-----------------------------------------------------------------------------
1432       subroutine make_var(n,idum,iter_csa)
1433
1434       use MD_data
1435       use energy_data
1436       use compare_data
1437       use control_data, only: vdisulf
1438       use geometry_data
1439       use geometry, only: dist
1440       include 'mpif.h'
1441
1442 !      implicit real*8 (a-h,o-z)
1443 !      include 'DIMENSIONS'
1444 !      include 'COMMON.IOUNITS'
1445 !      include 'COMMON.CSA'
1446 !      include 'COMMON.BANK'
1447 !      include 'COMMON.CHAIN'
1448 !      include 'COMMON.INTERACT'
1449 !      include 'COMMON.HAIRPIN'
1450 !      include 'COMMON.VAR'
1451 !      include 'COMMON.DISTFIT'
1452 !      include 'COMMON.GEO'
1453 !      include 'COMMON.CONTROL'
1454       logical :: nicht_getan,nicht_getan1,fail,lfound
1455       integer :: nharp,iharp(4,nres/3),nconf_harp
1456       integer :: iisucc(mxio)
1457       logical :: ifused(mxio)
1458       integer :: nhx_seed(nseed),ihx_seed(4,nres/3,nseed) !max_seed
1459       integer :: nhx_use(nseed),ihx_use(0:4,nres/3,nseed)
1460       integer :: nlx_seed(nseed),ilx_seed(2,nres/3,nseed),&
1461                nlx_use(nseed),ilx_use(nres/3,nseed)
1462 !      real(kind=4) :: ran1,ran2
1463
1464       integer :: i,j,k,n,idum,iter_csa,iran,index,n7frag,n8frag,n14frag,&
1465             n15frag,nbefrag,nlx_tot,iters,i1,i2,i3,ntot_gen,ngen,iih,&
1466             ij,jr,iim,nhx_tot,idummy,iter,iif,iig,icheck,ishift,iang,&
1467             n8c,ih_start,ih_end,n7c,index2,isize,nsucc,nacc,j1,nran,&
1468             ierror,ierrcode
1469       real(kind=8) :: d
1470
1471       write (iout,*) 'make_var : nseed=',nseed,'ntry=',n
1472       index=0
1473
1474 !-----------------------------------------
1475       if (n7.gt.0.or.n8.gt.0.or.n9.gt.0.or.n14.gt.0.or.n15.gt.0 &
1476           .or.n16.gt.0.or.n17.gt.0.or.n18.gt.0) &
1477                  call select_frag(n7frag,n8frag,n14frag,&
1478                  n15frag,nbefrag,iter_csa) 
1479
1480 !---------------------------------------------------
1481 ! N18 - random perturbation of one phi(=gamma) angle in a loop
1482 !
1483       IF (n18.gt.0) THEN
1484       nlx_tot=0
1485       do iters=1,nseed
1486         i1=is(iters)
1487         nlx_seed(iters)=0
1488         do i2=1,n14frag
1489           if (lvar_frag(i2,1).eq.i1) then
1490             nlx_seed(iters)=nlx_seed(iters)+5
1491             ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2)
1492             ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3)
1493             ilx_use(nlx_seed(iters),iters)=5
1494           endif
1495         enddo
1496         nlx_use(iters)=nlx_seed(iters)
1497         nlx_tot=nlx_tot+nlx_seed(iters)
1498       enddo
1499
1500       if (nlx_tot .ge. n18*nseed) then
1501         ntot_gen=n18*nseed
1502       else
1503         ntot_gen=(nlx_tot/nseed)*nseed
1504       endif
1505
1506       ngen=0
1507       do while (ngen.lt.ntot_gen)
1508       do iters=1,nseed
1509         iseed=is(iters)
1510         if (nlx_use(iters).gt.0) then
1511           nicht_getan=.true.
1512           do while (nicht_getan)
1513             iih=iran_num(1,nlx_seed(iters))
1514             if (ilx_use(iih,iters).gt.0) then
1515               nicht_getan=.false.
1516               ilx_use(iih,iters)=ilx_use(iih,iters)-1
1517               nlx_use(iters)=nlx_use(iters)-1
1518             endif
1519           enddo
1520           ngen=ngen+1
1521           index=index+1
1522           movenx(index)=18
1523           parent(1,index)=iseed
1524           parent(2,index)=0
1525
1526
1527            if (vdisulf) then
1528              nss_in(index)=bvar_nss(iseed)
1529              do ij=1,nss_in(index)
1530                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1531                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1532              enddo
1533            endif
1534            
1535
1536           do k=1,numch
1537             do j=2,nres-1
1538               do i=1,4
1539                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
1540               enddo
1541             enddo
1542           enddo
1543
1544           jr=iran_num(ilx_seed(1,iih,iters),ilx_seed(2,iih,iters))
1545           d=ran_number(-pi,pi)
1546           dihang_in(2,jr-2,1,index)=pinorm(dihang_in(2,jr-2,1,index)+d)
1547
1548
1549           if (ngen.eq.ntot_gen) goto 145
1550         endif
1551       enddo
1552       enddo
1553   145 continue
1554
1555       ENDIF
1556
1557
1558 !-----------------------------------------
1559 ! N17 : zip a beta in a seed by forcing one additional p-p contact
1560 !
1561       IF (n17.gt.0) THEN
1562       nhx_tot=0
1563       do iters=1,nseed
1564         i1=is(iters)
1565         nhx_seed(iters)=0
1566         nhx_use(iters)=0
1567         do i2=1,nbefrag
1568           if (avar_frag(i2,1).eq.i1) then
1569             nhx_seed(iters)=nhx_seed(iters)+1
1570             ihx_use(2,nhx_seed(iters),iters)=1
1571             if (avar_frag(i2,5)-avar_frag(i2,3).le.3.and. &
1572                  avar_frag(i2,2).gt.1.and.avar_frag(i2,4).lt.nres) then
1573              ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
1574              ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1
1575              ihx_use(0,nhx_seed(iters),iters)=1
1576              ihx_use(1,nhx_seed(iters),iters)=0
1577              nhx_use(iters)=nhx_use(iters)+1
1578             else
1579               if (avar_frag(i2,4).gt.avar_frag(i2,5)) then
1580                if (avar_frag(i2,2).gt.1.and. &
1581                            avar_frag(i2,4).lt.nres) then
1582                 ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
1583                 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1
1584                 ihx_use(0,nhx_seed(iters),iters)=1
1585                 ihx_use(1,nhx_seed(iters),iters)=0
1586                 nhx_use(iters)=nhx_use(iters)+1
1587                endif
1588                if (avar_frag(i2,3).lt.nres.and. &
1589                            avar_frag(i2,5).gt.1) then
1590                 ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1
1591                 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)-1
1592                 ihx_use(0,nhx_seed(iters),iters)= &
1593                           ihx_use(0,nhx_seed(iters),iters)+1
1594                 ihx_use(2,nhx_seed(iters),iters)=0
1595                 nhx_use(iters)=nhx_use(iters)+1
1596                endif
1597               else
1598                if (avar_frag(i2,2).gt.1.and. &
1599                            avar_frag(i2,4).gt.1) then
1600                 ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1
1601                 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)-1
1602                 ihx_use(0,nhx_seed(iters),iters)=1
1603                 ihx_use(1,nhx_seed(iters),iters)=0
1604                 nhx_use(iters)=nhx_use(iters)+1
1605                endif
1606                if (avar_frag(i2,3).lt.nres.and. &
1607                            avar_frag(i2,5).lt.nres) then
1608                 ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1
1609                 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)+1
1610                 ihx_use(0,nhx_seed(iters),iters)= &
1611                           ihx_use(0,nhx_seed(iters),iters)+1
1612                 ihx_use(2,nhx_seed(iters),iters)=0
1613                 nhx_use(iters)=nhx_use(iters)+1
1614                endif
1615               endif
1616             endif
1617           endif
1618         enddo
1619
1620         nhx_tot=nhx_tot+nhx_use(iters)
1621 !d        write (iout,*) "debug N17",iters,nhx_seed(iters),
1622 !d     &                     nhx_use(iters),nhx_tot
1623       enddo
1624
1625       if (nhx_tot .ge. n17*nseed) then
1626         ntot_gen=n17*nseed
1627       else if (nhx_tot .ge. nseed) then
1628         ntot_gen=(nhx_tot/nseed)*nseed
1629       else
1630         ntot_gen=nhx_tot
1631       endif
1632 !d      write (iout,*) "debug N17==",ntot_gen,nhx_tot,nseed
1633
1634       ngen=0
1635       do while (ngen.lt.ntot_gen)
1636       do iters=1,nseed
1637         iseed=is(iters)
1638         if (nhx_use(iters).gt.0) then
1639 !d          write (iout,*) "debug N17",nhx_use(iters),ngen,ntot_gen
1640 !d          write (iout,*) "debugN17^",
1641 !d     &                    (ihx_use(0,k,iters),k=1,nhx_use(iters))
1642           nicht_getan=.true.
1643           do while (nicht_getan)
1644             iih=iran_num(1,nhx_seed(iters))
1645 !d            write (iout,*) "debugN17^",iih
1646             if (ihx_use(0,iih,iters).gt.0) then
1647               iim=iran_num(1,2)
1648 !d                write (iout,*) "debugN17=",iih,nhx_seed(iters)
1649 !d                write (iout,*) "debugN17-",iim,'##',
1650 !d     &                           (ihx_use(k,iih,iters),k=0,2)
1651 !d                call flush(iout)
1652               do while (ihx_use(iim,iih,iters).eq.1)
1653                 iim=iran_num(1,2)
1654 !d                write (iout,*) "debugN17-",iim,'##',
1655 !d     &                           (ihx_use(k,iih,iters),k=0,2)
1656 !d                call flush(iout)
1657               enddo
1658               nicht_getan=.false.
1659               ihx_use(iim,iih,iters)=1
1660               ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 
1661               nhx_use(iters)=nhx_use(iters)-1
1662             endif
1663           enddo
1664           ngen=ngen+1
1665           index=index+1
1666           movenx(index)=17
1667           parent(1,index)=iseed
1668           parent(2,index)=0
1669
1670            if (vdisulf) then
1671              nss_in(index)=bvar_nss(iseed)
1672              do ij=1,nss_in(index)
1673                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1674                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1675              enddo
1676            endif
1677
1678           do k=1,numch
1679             do j=2,nres-1
1680               do i=1,4
1681                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
1682               enddo
1683             enddo
1684           enddo
1685
1686           if (iim.eq.1) then
1687             idata(1,index)=ihx_seed(1,iih,iters)
1688             idata(2,index)=ihx_seed(2,iih,iters)
1689           else
1690             idata(1,index)=ihx_seed(3,iih,iters)
1691             idata(2,index)=ihx_seed(4,iih,iters)
1692           endif
1693
1694           if (ngen.eq.ntot_gen) goto 115
1695         endif
1696       enddo
1697       enddo
1698   115 continue
1699       write (iout,*) "N17",n17," ngen/nseed",ngen/nseed,&
1700                                  ngen,nseed
1701
1702
1703       ENDIF
1704 !-----------------------------------------
1705 ! N16 : slide non local beta in a seed by +/- 1 or +/- 2
1706 !
1707       IF (n16.gt.0) THEN
1708       nhx_tot=0
1709       do iters=1,nseed
1710         i1=is(iters)
1711         nhx_seed(iters)=0
1712         do i2=1,n7frag
1713           if (bvar_frag(i2,1).eq.i1) then
1714             nhx_seed(iters)=nhx_seed(iters)+1
1715             ihx_seed(1,nhx_seed(iters),iters)=bvar_frag(i2,3)
1716             ihx_seed(2,nhx_seed(iters),iters)=bvar_frag(i2,4)
1717             ihx_seed(3,nhx_seed(iters),iters)=bvar_frag(i2,5)
1718             ihx_seed(4,nhx_seed(iters),iters)=bvar_frag(i2,6)
1719             ihx_use(0,nhx_seed(iters),iters)=4
1720             do i3=1,4
1721               ihx_use(i3,nhx_seed(iters),iters)=0
1722             enddo
1723           endif
1724         enddo
1725         nhx_use(iters)=4*nhx_seed(iters)
1726         nhx_tot=nhx_tot+nhx_seed(iters)
1727 !d        write (iout,*) "debug N16",iters,nhx_seed(iters)
1728       enddo
1729
1730       if (4*nhx_tot .ge. n16*nseed) then
1731         ntot_gen=n16*nseed
1732       else if (4*nhx_tot .ge. nseed) then
1733         ntot_gen=(4*nhx_tot/nseed)*nseed
1734       else
1735         ntot_gen=4*nhx_tot
1736       endif
1737       write (iout,*) "debug N16",ntot_gen,4*nhx_tot,nseed
1738
1739       ngen=0
1740       do while (ngen.lt.ntot_gen)
1741       do iters=1,nseed
1742         iseed=is(iters)
1743         if (nhx_use(iters).gt.0) then
1744           nicht_getan=.true.
1745           do while (nicht_getan)
1746             iih=iran_num(1,nhx_seed(iters))
1747             if (ihx_use(0,iih,iters).gt.0) then
1748               iim=iran_num(1,4)
1749               do while (ihx_use(iim,iih,iters).eq.1)
1750 !d                write (iout,*) iim,
1751 !d     &                ihx_use(0,iih,iters),ihx_use(iim,iih,iters)
1752                 iim=iran_num(1,4)
1753               enddo
1754               nicht_getan=.false.
1755               ihx_use(iim,iih,iters)=1
1756               ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 
1757               nhx_use(iters)=nhx_use(iters)-1
1758             endif
1759           enddo
1760           ngen=ngen+1
1761           index=index+1
1762           movenx(index)=16
1763           parent(1,index)=iseed
1764           parent(2,index)=0
1765
1766            if (vdisulf) then
1767              nss_in(index)=bvar_nss(iseed)
1768              do ij=1,nss_in(index)
1769                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1770                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1771              enddo
1772            endif
1773
1774           do k=1,numch
1775             do j=2,nres-1
1776               do i=1,4
1777                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
1778               enddo
1779             enddo
1780           enddo
1781
1782           do i=1,4
1783            idata(i,index)=ihx_seed(i,iih,iters)
1784           enddo
1785           idata(5,index)=iim           
1786
1787           if (ngen.eq.ntot_gen) goto 116
1788         endif
1789       enddo
1790       enddo
1791   116 continue
1792       write (iout,*) "N16",n16," ngen/nseed",ngen/nseed,&
1793                                  ngen,nseed
1794       ENDIF
1795 !-----------------------------------------
1796 ! N15 : copy two 2nd structure elements from 1 or 2 conf. in bank to a seed
1797 !
1798       IF (n15.gt.0) THEN
1799
1800        do iters=1,nseed
1801          iseed=is(iters)
1802          do i=1,mxio
1803            ifused(i)=.false.
1804          enddo
1805
1806          do idummy=1,n15
1807            iter=0
1808    84      continue
1809
1810            iran=0
1811            iif=iran_num(1,n15frag)
1812            do while( (ifused(iif) .or. svar_frag(iif,1).eq.iseed) .and. &
1813                             iran.le.mxio )
1814             iif=iran_num(1,n15frag)
1815             iran=iran+1
1816            enddo
1817            if(iran.ge.mxio) goto 811
1818
1819            iran=0
1820            iig=iran_num(1,n15frag)
1821            do while( (ifused(iig) .or. svar_frag(iig,1).eq.iseed  .or. &
1822                    .not.(svar_frag(iif,3).lt.svar_frag(iig,2).or. &
1823                          svar_frag(iig,3).lt.svar_frag(iif,2)) ) .and. &
1824                             iran.le.mxio )
1825             iig=iran_num(1,n15frag)
1826             iran=iran+1
1827            enddo
1828            if(iran.ge.mxio) goto 811
1829
1830            index=index+1
1831            movenx(index)=15
1832            parent(1,index)=iseed
1833            parent(2,index)=svar_frag(iif,1)
1834            parent(3,index)=svar_frag(iig,1)
1835
1836
1837            if (vdisulf) then
1838              nss_in(index)=bvar_nss(iseed)
1839              do ij=1,nss_in(index)
1840                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1841                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1842              enddo
1843            endif
1844
1845            ifused(iif)=.true.
1846            ifused(iig)=.true.
1847            call newconf_copy(idum,dihang_in(1,1,1,index),&
1848                    svar_frag(iif,1),svar_frag(iif,2),svar_frag(iif,3))
1849
1850            do j=svar_frag(iig,2),svar_frag(iig,3)
1851              do i=1,4
1852               dihang_in(i,j,1,index)=bvar(i,j,1,svar_frag(iig,1))
1853              enddo
1854            enddo
1855
1856
1857            if(iter.lt.10) then
1858             call check_old(icheck,index)
1859             if(icheck.eq.1) then 
1860               index=index-1
1861               ifused(iif)=.false. 
1862               goto 84
1863             endif
1864            endif
1865
1866   811     continue
1867          enddo
1868        enddo
1869       ENDIF
1870
1871 !-----------------------------------------
1872 ! N14 local_move (Maurizio) for loops in a seed
1873 !
1874       IF (n14.gt.0) THEN
1875       nlx_tot=0
1876       do iters=1,nseed
1877         i1=is(iters)
1878         nlx_seed(iters)=0
1879         do i2=1,n14frag
1880           if (lvar_frag(i2,1).eq.i1) then
1881             nlx_seed(iters)=nlx_seed(iters)+3
1882             ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2)
1883             ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3)
1884             ilx_use(nlx_seed(iters),iters)=3
1885           endif
1886         enddo
1887         nlx_use(iters)=nlx_seed(iters)
1888         nlx_tot=nlx_tot+nlx_seed(iters)
1889 !d        write (iout,*) "debug N14",iters,nlx_seed(iters)
1890       enddo
1891
1892       if (nlx_tot .ge. n14*nseed) then
1893         ntot_gen=n14*nseed
1894       else
1895         ntot_gen=(nlx_tot/nseed)*nseed
1896       endif
1897 !d      write (iout,*) "debug N14",ntot_gen,n14frag,nseed
1898
1899       ngen=0
1900       do while (ngen.lt.ntot_gen)
1901       do iters=1,nseed
1902         iseed=is(iters)
1903         if (nlx_use(iters).gt.0) then
1904           nicht_getan=.true.
1905           do while (nicht_getan)
1906             iih=iran_num(1,nlx_seed(iters))
1907             if (ilx_use(iih,iters).gt.0) then
1908               nicht_getan=.false.
1909               ilx_use(iih,iters)=ilx_use(iih,iters)-1
1910               nlx_use(iters)=nlx_use(iters)-1
1911             endif
1912           enddo
1913           ngen=ngen+1
1914           index=index+1
1915           movenx(index)=14
1916           parent(1,index)=iseed
1917           parent(2,index)=0
1918
1919           idata(1,index)=ilx_seed(1,iih,iters)
1920           idata(2,index)=ilx_seed(2,iih,iters)
1921
1922
1923            if (vdisulf) then
1924              nss_in(index)=bvar_nss(iseed)
1925              do ij=1,nss_in(index)
1926                iss_in(ij,index)=bvar_ss(1,ij,iseed)
1927                jss_in(ij,index)=bvar_ss(2,ij,iseed)
1928              enddo
1929            endif
1930            
1931
1932           do k=1,numch
1933             do j=2,nres-1
1934               do i=1,4
1935                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
1936               enddo
1937             enddo
1938           enddo
1939
1940           if (ngen.eq.ntot_gen) goto 131
1941         endif
1942       enddo
1943       enddo
1944   131 continue
1945 !d      write (iout,*) "N14",n14," ngen/nseed",ngen/nseed,
1946 !d     &                           ngen,nseed
1947
1948       ENDIF
1949 !-----------------------------------------
1950 ! N9 : shift a helix in a seed
1951 !
1952       IF (n9.gt.0) THEN
1953       nhx_tot=0
1954       do iters=1,nseed
1955         i1=is(iters)
1956         nhx_seed(iters)=0
1957         do i2=1,n8frag
1958           if (hvar_frag(i2,1).eq.i1) then
1959             nhx_seed(iters)=nhx_seed(iters)+1
1960             ihx_seed(1,nhx_seed(iters),iters)=hvar_frag(i2,2)
1961             ihx_seed(2,nhx_seed(iters),iters)=hvar_frag(i2,3)
1962             ihx_use(0,nhx_seed(iters),iters)=4
1963             do i3=1,4
1964               ihx_use(i3,nhx_seed(iters),iters)=0
1965             enddo
1966           endif
1967         enddo
1968         nhx_use(iters)=4*nhx_seed(iters)
1969         nhx_tot=nhx_tot+nhx_seed(iters)
1970 !d        write (iout,*) "debug N9",iters,nhx_seed(iters)
1971       enddo
1972
1973       if (4*nhx_tot .ge. n9*nseed) then
1974         ntot_gen=n9*nseed
1975       else
1976         ntot_gen=(4*nhx_tot/nseed)*nseed
1977       endif
1978 !d      write (iout,*) "debug N9",ntot_gen,n8frag,nseed
1979
1980       ngen=0
1981       do while (ngen.lt.ntot_gen)
1982       do iters=1,nseed
1983         iseed=is(iters)
1984         if (nhx_use(iters).gt.0) then
1985           nicht_getan=.true.
1986           do while (nicht_getan)
1987             iih=iran_num(1,nhx_seed(iters))
1988             if (ihx_use(0,iih,iters).gt.0) then
1989               iim=iran_num(1,4)
1990               do while (ihx_use(iim,iih,iters).eq.1)
1991 !d                write (iout,*) iim,
1992 !d     &                ihx_use(0,iih,iters),ihx_use(iim,iih,iters)
1993                 iim=iran_num(1,4)
1994               enddo
1995               nicht_getan=.false.
1996               ihx_use(iim,iih,iters)=1
1997               ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 
1998               nhx_use(iters)=nhx_use(iters)-1
1999             endif
2000           enddo
2001           ngen=ngen+1
2002           index=index+1
2003           movenx(index)=9
2004           parent(1,index)=iseed
2005           parent(2,index)=0
2006
2007            if (vdisulf) then
2008              nss_in(index)=bvar_nss(iseed)
2009              do ij=1,nss_in(index)
2010                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2011                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2012              enddo
2013            endif
2014
2015           do k=1,numch
2016             do j=2,nres-1
2017               do i=1,4
2018                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2019               enddo
2020             enddo
2021           enddo
2022
2023           jstart=max(nnt,ihx_seed(1,iih,iters)+1)
2024           jend=min(nct,ihx_seed(2,iih,iters))
2025 !d          write (iout,*) "debug N9",iters,iih,jstart,jend
2026           if (iim.eq.1) then
2027             ishift=-2
2028           else if (iim.eq.2) then
2029             ishift=-1
2030           else if (iim.eq.3) then
2031             ishift=1
2032           else if (iim.eq.4) then
2033             ishift=2
2034           else
2035             write (iout,*) 'CHUJ NASTAPIL: iim=',iim
2036 #ifdef MPI !el
2037             call mpi_abort(mpi_comm_world,ierror,ierrcode)
2038 #endif
2039           endif
2040           do j=jstart,jend
2041             if (itype(j).eq.10) then
2042               iang=2
2043             else
2044               iang=4
2045             endif
2046             do i=1,iang
2047               if (j+ishift.ge.nnt.and.j+ishift.le.nct) &
2048                 dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed)
2049             enddo
2050           enddo
2051           if (ishift.gt.0) then
2052            do j=0,ishift-1
2053             if (itype(jend+j).eq.10) then
2054               iang=2
2055             else
2056               iang=4
2057             endif
2058             do i=1,iang
2059               if (jend+j.ge.nnt.and.jend+j.le.nct) &
2060                 dihang_in(i,jstart+j,1,index)=bvar(i,jend+j,1,iseed)
2061             enddo
2062            enddo
2063           else
2064            do j=0,-ishift-1
2065             if (itype(jstart+j).eq.10) then
2066               iang=2
2067             else
2068               iang=4
2069             endif
2070             do i=1,iang
2071               if (jend+j.ge.nnt.and.jend+j.le.nct) &
2072                 dihang_in(i,jend+j,1,index)=bvar(i,jstart+j,1,iseed)
2073             enddo
2074            enddo
2075           endif
2076           if (ngen.eq.ntot_gen) goto 133
2077         endif
2078       enddo
2079       enddo
2080   133 continue
2081 !d      write (iout,*) "N9",n9," ngen/nseed",ngen/nseed,
2082 !d     &                           ngen,nseed
2083
2084       ENDIF
2085 !-----------------------------------------
2086 ! N8 : copy a helix from bank to seed
2087 !
2088       if (n8.gt.0) then 
2089        if (n8frag.lt.n8) then 
2090           write (iout,*) "N8: only ",n8frag,'helices'
2091           n8c=n8frag
2092        else
2093           n8c=n8
2094        endif
2095
2096        do iters=1,nseed
2097          iseed=is(iters)
2098          do i=1,mxio
2099            ifused(i)=.false.
2100          enddo
2101
2102
2103          do idummy=1,n8c
2104            iter=0
2105    94      continue
2106            iran=0
2107            iif=iran_num(1,n8frag)
2108            do while( (ifused(iif) .or. hvar_frag(iif,1).eq.iseed) .and. &
2109                             iran.le.mxio )
2110             iif=iran_num(1,n8frag)
2111             iran=iran+1
2112            enddo
2113             
2114            if(iran.ge.mxio) goto 911
2115
2116            index=index+1
2117            movenx(index)=8
2118            parent(1,index)=iseed
2119            parent(2,index)=hvar_frag(iif,1)
2120
2121
2122            if (vdisulf) then
2123              nss_in(index)=bvar_nss(iseed)
2124              do ij=1,nss_in(index)
2125                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2126                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2127              enddo
2128            endif
2129
2130            ifused(iif)=.true.
2131            if (hvar_frag(iif,3)-hvar_frag(iif,2).le.6) then
2132              call newconf_copy(idum,dihang_in(1,1,1,index),&
2133                    hvar_frag(iif,1),hvar_frag(iif,2),hvar_frag(iif,3))
2134            else
2135               ih_start=iran_num(hvar_frag(iif,2),hvar_frag(iif,3)-6)
2136               ih_end=iran_num(ih_start,hvar_frag(iif,3))
2137              call newconf_copy(idum,dihang_in(1,1,1,index),&
2138                    hvar_frag(iif,1),ih_start,ih_end)
2139            endif
2140            iter=iter+1
2141            if(iter.lt.10) then
2142             call check_old(icheck,index)
2143             if(icheck.eq.1) then 
2144               index=index-1
2145               ifused(iif)=.false. 
2146               goto 94
2147             endif
2148            endif
2149
2150
2151   911     continue
2152
2153          enddo
2154        enddo
2155
2156       endif
2157
2158 !-----------------------------------------
2159 ! N7 : copy nonlocal beta fragment from bank to seed
2160 !
2161       if (n7.gt.0) then 
2162        if (n7frag.lt.n7) then 
2163           write (iout,*) "N7: only ",n7frag,'nonlocal fragments'
2164           n7c=n7frag
2165        else
2166           n7c=n7
2167        endif
2168
2169        do i=1,nres
2170          do j=1,mxio2
2171             iff_in(i,j)=0
2172          enddo
2173        enddo
2174        index2=0
2175        do i=1,mxio
2176          isend2(i)=0
2177        enddo
2178
2179        do iters=1,nseed
2180          iseed=is(iters)
2181          do i=1,mxio
2182            ifused(i)=.false.
2183          enddo
2184
2185          do idummy=1,n7c
2186            iran=0
2187            iif=iran_num(1,n7frag)
2188            do while( (ifused(iif) .or. bvar_frag(iif,1).eq.iseed) .and. &
2189                             iran.le.mxio )
2190             iif=iran_num(1,n7frag)
2191             iran=iran+1
2192            enddo
2193             
2194 !d       write (*,'(3i5,l,4i5)'),iters,idummy,iif,ifused(iif),
2195 !d     &                    bvar_frag(iif,1),iseed,iran,index2
2196          
2197            if(iran.ge.mxio) goto 999
2198            if(index2.ge.mxio2) goto 999 
2199
2200            index=index+1
2201            movenx(index)=7
2202            parent(1,index)=iseed
2203            parent(2,index)=bvar_frag(iif,1)
2204            index2=index2+1
2205            isend2(index)=index2
2206            ifused(iif)=.true.
2207
2208            if (vdisulf) then
2209              nss_in(index)=bvar_nss(iseed)
2210              do ij=1,nss_in(index)
2211                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2212                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2213              enddo
2214            endif
2215
2216            do k=1,numch
2217             do j=2,nres-1
2218              do i=1,4
2219               dihang_in2(i,j,k,index2)=bvar(i,j,k,bvar_frag(iif,1))
2220              enddo
2221             enddo
2222            enddo
2223
2224            if (bvar_frag(iif,2).eq.4) then                                                                  
2225              do i=bvar_frag(iif,3),bvar_frag(iif,4)
2226                iff_in(i,index2)=1                                                              
2227              enddo                                                                   
2228              if (bvar_frag(iif,5).lt.bvar_frag(iif,6)) then
2229 !d               print *,'###',bvar_frag(iif,3),bvar_frag(iif,4),
2230 !d     &                 bvar_frag(iif,5),bvar_frag(iif,6)
2231                do i=bvar_frag(iif,5),bvar_frag(iif,6)
2232                  iff_in(i,index2)=1                                                              
2233                enddo                                                                   
2234              else
2235 !d               print *,'###',bvar_frag(iif,3),bvar_frag(iif,4),
2236 !d     &                 bvar_frag(iif,6),bvar_frag(iif,5)
2237                do i=bvar_frag(iif,6),bvar_frag(iif,5)
2238                  iff_in(i,index2)=1                                                              
2239                enddo                                                                   
2240              endif
2241            endif
2242
2243            do k=1,numch
2244             do j=2,nres-1
2245              do i=1,4
2246               dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2247              enddo
2248             enddo
2249            enddo
2250
2251
2252   999     continue
2253
2254          enddo
2255        enddo
2256
2257       endif
2258 !-----------------------------------------------
2259 ! N6 : copy random continues fragment from bank to seed
2260 !
2261       do iters=1,nseed
2262        iseed=is(iters)
2263        do idummy=1,n6
2264         isize=(is2-is1+1)*ran1(idum)+is1
2265         index=index+1
2266         movenx(index)=6
2267
2268
2269            if (vdisulf) then
2270              nss_in(index)=bvar_nss(iseed)
2271              do ij=1,nss_in(index)
2272                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2273                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2274              enddo
2275            endif
2276
2277         iter=0
2278   104   continue
2279         if(icycle.le.0) then
2280          i1=nconf*  ran1(idum)+1
2281          i1=nbank-nconf+i1
2282         else
2283          i1=nbank*  ran1(idum)+1
2284         endif
2285         if(i1.eq.iseed) goto 104
2286         iter=iter+1
2287         call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
2288         parent(1,index)=iseed
2289         parent(2,index)=i1
2290         if(iter.lt.10) then
2291          call check_old(icheck,index)
2292          if(icheck.eq.1) goto 104
2293         endif
2294        enddo
2295       enddo
2296 !-----------------------------------------
2297       if (n3.gt.0.or.n4.gt.0) call gen_hairpin
2298       nconf_harp=0
2299       do iters=1,nseed
2300         if (nharp_seed(iters).gt.0) nconf_harp=nconf_harp+1
2301       enddo      
2302 !-----------------------------------------
2303 ! N3 : copy hairpin from bank to seed
2304 !
2305       do iters=1,nseed
2306        iseed=is(iters)
2307        nsucc=0
2308        nacc=0
2309        do idummy=1,n3
2310         index=index+1
2311         iter=0
2312   124   continue
2313         if(icycle.le.0) then
2314          i1=nconf*  ran1(idum)+1
2315          i1=nbank-nconf+i1
2316         else
2317          i1=nbank*  ran1(idum)+1
2318         endif
2319         if(i1.eq.iseed) goto 124
2320         do k=1,nsucc
2321           if (i1.eq.iisucc(k).and.nsucc.lt.nconf_harp-1) goto 124
2322         enddo
2323         nsucc=nsucc+1
2324         iisucc(nsucc)=i1
2325         iter=iter+1
2326         call newconf_residue_hairpin(idum,dihang_in(1,1,1,index),&
2327            i1,fail)
2328         if (fail) then
2329           if (icycle.le.0 .and. nsucc.eq.nconf .or. &
2330               icycle.gt.0 .and. nsucc.eq.nbank)  then
2331             index=index-1
2332             goto 125
2333           else
2334             goto 124
2335           endif
2336         endif
2337         if(iter.lt.10) then
2338          call check_old(icheck,index)
2339          if(icheck.eq.1) goto 124
2340         endif
2341         movenx(index)=3
2342         parent(1,index)=iseed
2343         parent(2,index)=i1
2344
2345
2346            if (vdisulf) then
2347              nss_in(index)=bvar_nss(iseed)
2348              do ij=1,nss_in(index)
2349                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2350                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2351              enddo
2352            endif
2353
2354         nacc=nacc+1
2355        enddo
2356 ! if not enough hairpins, supplement with windows
2357   125  continue
2358 !dd       if (n3.ne.0) write (iout,*) "N3",n3," nsucc",nsucc," nacc",nacc 
2359        do idummy=nacc+1,n3
2360         isize=(is2-is1+1)*ran1(idum)+is1
2361         index=index+1
2362         movenx(index)=6
2363         parent(1,index)=iseed
2364         parent(2,index)=i1
2365
2366
2367            if (vdisulf) then
2368              nss_in(index)=bvar_nss(iseed)
2369              do ij=1,nss_in(index)
2370                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2371                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2372              enddo
2373            endif
2374
2375         iter=0
2376   114   continue
2377         if(icycle.le.0) then
2378          i1=nconf*  ran1(idum)+1
2379          i1=nbank-nconf+i1
2380         else
2381          i1=nbank*  ran1(idum)+1
2382         endif
2383         if(i1.eq.iseed) goto 114
2384         iter=iter+1
2385         call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
2386         if(iter.lt.10) then
2387          call check_old(icheck,index)
2388          if(icheck.eq.1) goto 114
2389         endif
2390        enddo
2391       enddo
2392 !-----------------------------------------
2393 ! N4 : shift a turn in hairpin in seed
2394 !
2395       IF (N4.GT.0) THEN
2396       if (4*nharp_tot .ge. n4*nseed) then
2397         ntot_gen=n4*nseed
2398       else
2399         ntot_gen=(4*nharp_tot/nseed)*nseed
2400       endif
2401       ngen=0
2402       do while (ngen.lt.ntot_gen)
2403       do iters=1,nseed
2404         iseed=is(iters)
2405 !        write (iout,*) 'iters',iters,' iseed',iseed,' nharp_seed',
2406 !     &     nharp_seed(iters),' nharp_use',nharp_use(iters),
2407 !     &     ' ntot_gen',ntot_gen
2408 !        write (iout,*) 'iharp_use(0)',
2409 !     &     (iharp_use(0,k,iters),k=1,nharp_seed(iters))
2410         if (nharp_use(iters).gt.0) then
2411           nicht_getan=.true.
2412           do while (nicht_getan)
2413             iih=iran_num(1,nharp_seed(iters))
2414 !            write (iout,*) 'iih',iih,' iharp_use',
2415 !     &            (iharp_use(k,iih,iters),k=1,4)
2416             if (iharp_use(0,iih,iters).gt.0) then
2417               nicht_getan1=.true.
2418               do while (nicht_getan1)
2419                 iim=iran_num(1,4)
2420                 nicht_getan1=iharp_use(iim,iih,iters).eq.1
2421               enddo
2422               nicht_getan=.false.
2423               iharp_use(iim,iih,iters)=1
2424               iharp_use(0,iih,iters)=iharp_use(0,iih,iters)-1 
2425               nharp_use(iters)=nharp_use(iters)-1
2426 !dd              write (iout,'(a16,i3,a5,i2,a10,2i4)') 
2427 !dd     &           'N4 selected hairpin',iih,' move',iim,' iharp_seed',
2428 !dd     &            iharp_seed(1,iih,iters),iharp_seed(2,iih,iters)
2429             endif
2430           enddo
2431           ngen=ngen+1
2432           index=index+1
2433           movenx(index)=4
2434           parent(1,index)=iseed
2435           parent(2,index)=0
2436
2437
2438            if (vdisulf) then
2439              nss_in(index)=bvar_nss(iseed)
2440              do ij=1,nss_in(index)
2441                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2442                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2443              enddo
2444            endif
2445
2446           do k=1,numch
2447             do j=2,nres-1
2448               do i=1,4
2449                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2450               enddo
2451             enddo
2452           enddo
2453           jstart=iharp_seed(1,iih,iters)+1
2454           jend=iharp_seed(2,iih,iters)
2455           if (iim.eq.1) then
2456             ishift=-2
2457           else if (iim.eq.2) then
2458             ishift=-1
2459           else if (iim.eq.3) then
2460             ishift=1
2461           else if (iim.eq.4) then
2462             ishift=2
2463           else
2464             write (iout,*) 'CHUJ NASTAPIL: iim=',iim
2465 #ifdef MPI !el
2466             call mpi_abort(mpi_comm_world,ierror,ierrcode)
2467 #endif !el
2468           endif
2469 !          write (iout,*) 'jstart',jstart,' jend',jend,' ishift',ishift
2470 !          write (iout,*) 'Before turn shift'
2471 !          do j=2,nres-1
2472 !            theta(j+1)=dihang_in(1,j,1,index)
2473 !            phi(j+2)=dihang_in(2,j,1,index)
2474 !            alph(j)=dihang_in(3,j,1,index)
2475 !            omeg(j)=dihang_in(4,j,1,index)
2476 !          enddo
2477 !          call intout
2478           do j=jstart,jend
2479             if (itype(j).eq.10) then
2480               iang=2
2481             else
2482               iang=4
2483             endif
2484             do i=1,iang
2485               if (j+ishift.ge.nnt.and.j+ishift.le.nct) &
2486                 dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed)
2487             enddo
2488           enddo
2489 !          write (iout,*) 'After turn shift'
2490 !          do j=2,nres-1
2491 !            theta(j+1)=dihang_in(1,j,1,index)
2492 !            phi(j+2)=dihang_in(2,j,1,index)
2493 !            alph(j)=dihang_in(3,j,1,index)
2494 !            omeg(j)=dihang_in(4,j,1,index)
2495 !          enddo
2496 !          call intout
2497           if (ngen.eq.ntot_gen) goto 135
2498         endif
2499       enddo
2500       enddo
2501 ! if not enough hairpins, supplement with windows
2502 !      write (iout,*) 'end of enddo'
2503   135 continue
2504 !dd      write (iout,*) "N4",n4," ngen/nseed",ngen/nseed,
2505 !dd     &                           ngen,nseed
2506       do iters=1,nseed
2507         iseed=is(iters)
2508         do idummy=ngen/nseed+1,n4
2509           isize=(is2-is1+1)*ran1(idum)+is1
2510           index=index+1
2511           movenx(index)=6
2512
2513            if (vdisulf) then
2514              nss_in(index)=bvar_nss(iseed)
2515              do ij=1,nss_in(index)
2516                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2517                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2518              enddo
2519            endif
2520
2521
2522           iter=0
2523   134     continue
2524           if(icycle.le.0) then
2525           i1=nconf*  ran1(idum)+1
2526           i1=nbank-nconf+i1
2527           else
2528             i1=nbank*  ran1(idum)+1
2529           endif
2530           if(i1.eq.iseed) goto 134
2531             iter=iter+1
2532             call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
2533             parent(1,index)=iseed
2534             parent(2,index)=i1
2535             if(iter.lt.10) then
2536               call check_old(icheck,index)
2537             if(icheck.eq.1) goto 134
2538           endif
2539         enddo
2540       enddo
2541       ENDIF
2542 !-----------------------------------------
2543 ! N5 : copy one residue from bank to seed (normally switched off - use N1)
2544 !
2545       do iters=1,nseed
2546        iseed=is(iters)
2547        isize=1
2548        do i=1,n5
2549         index=index+1
2550         movenx(index)=5
2551
2552            if (vdisulf) then
2553              nss_in(index)=bvar_nss(iseed)
2554              do ij=1,nss_in(index)
2555                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2556                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2557              enddo
2558            endif
2559
2560
2561         iter=0
2562   105   continue
2563         if(icycle.le.0) then
2564          i1=nconf*  ran1(idum)+1
2565          i1=nbank-nconf+i1
2566         else
2567          i1=nbank*  ran1(idum)+1
2568         endif
2569         if(i1.eq.iseed) goto 105
2570         iter=iter+1
2571         call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize)
2572         parent(1,index)=iseed
2573         parent(2,index)=i1
2574         if(iter.lt.10) then
2575          call check_old(icheck,index)
2576          if(icheck.eq.1) goto 105
2577         endif
2578        enddo
2579       enddo
2580 !-----------------------------------------
2581 ! N2 : copy backbone of one residue from bank or first bank to seed
2582 ! (normally switched off - use N1)
2583 !
2584       do iters=1,nseed
2585        iseed=is(iters)
2586        do i=n2,1,-1
2587         if(icycle.le.0.and.iuse.gt.nconf-irr) then
2588          iseed=ran1(idum)*nconf+1
2589          iseed=nbank-nconf+iseed
2590         endif
2591         index=index+1
2592         movenx(index)=2
2593
2594         if (vdisulf) then
2595              nss_in(index)=bvar_nss(iseed)
2596              do ij=1,nss_in(index)
2597                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2598                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2599              enddo
2600         endif
2601
2602         iter=0
2603   102   i1=  ran1(idum)*nbank+1
2604         if(i1.eq.iseed) goto 102
2605         iter=iter+1
2606         if(icycle.le.0.and.iuse.gt.nconf-irr) then
2607          nran=mod(i-1,nran0)+3
2608          call newconf1arr(idum,dihang_in(1,1,1,index),nran,i1)
2609          parent(1,index)=-iseed
2610          parent(2,index)=-i1
2611         else if(icycle.le.0.and.iters.le.iuse) then
2612          nran=mod(i-1,nran0)+1
2613          call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1)
2614          parent(1,index)=iseed
2615          parent(2,index)=-i1
2616         else
2617          nran=mod(i-1,nran1)+1
2618          if(ran1(idum).lt.0.5) then
2619           call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1)
2620           parent(1,index)=iseed
2621           parent(2,index)=-i1
2622          else
2623           call newconf1abb(idum,dihang_in(1,1,1,index),nran,i1)
2624           parent(1,index)=iseed
2625           parent(2,index)=i1
2626          endif
2627         endif
2628         if(iter.lt.10) then
2629          call check_old(icheck,index)
2630          if(icheck.eq.1) goto 102
2631         endif
2632        enddo
2633       enddo
2634 !-----------------------------------------
2635 ! N1 : copy backbone or sidechain of one residue from bank or 
2636 !      first bank to seed
2637 !
2638       do iters=1,nseed
2639        iseed=is(iters)
2640        do i=n1,1,-1
2641         if(icycle.le.0.and.iuse.gt.nconf-irr) then
2642          iseed=ran1(idum)*nconf+1
2643          iseed=nbank-nconf+iseed
2644         endif
2645         index=index+1
2646         movenx(index)=1
2647
2648         if (vdisulf) then
2649              nss_in(index)=bvar_nss(iseed)
2650              do ij=1,nss_in(index)
2651                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2652                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2653              enddo
2654         endif
2655
2656         iter=0
2657   101   i1=  ran1(idum)*nbank+1
2658
2659         if(i1.eq.iseed) goto 101
2660         iter=iter+1
2661         if(icycle.le.0.and.iuse.gt.nconf-irr) then
2662          nran=mod(i-1,nran0)+3
2663          call newconf1rr(idum,dihang_in(1,1,1,index),nran,i1)
2664          parent(1,index)=-iseed
2665          parent(2,index)=-i1
2666         else if(icycle.le.0.and.iters.le.iuse) then
2667          nran=mod(i-1,nran0)+1
2668          call newconf1br(idum,dihang_in(1,1,1,index),nran,i1)
2669          parent(1,index)=iseed
2670          parent(2,index)=-i1
2671         else
2672          nran=mod(i-1,nran1)+1
2673          if(ran1(idum).lt.0.5) then
2674           call newconf1br(idum,dihang_in(1,1,1,index),nran,i1)
2675           parent(1,index)=iseed
2676           parent(2,index)=-i1
2677          else
2678           call newconf1bb(idum,dihang_in(1,1,1,index),nran,i1)
2679           parent(1,index)=iseed
2680           parent(2,index)=i1
2681          endif
2682         endif
2683         if(iter.lt.10) then
2684          call check_old(icheck,index)
2685          if(icheck.eq.1) goto 101
2686         endif
2687        enddo
2688       enddo
2689 !-----------------------------------------
2690 ! N0 just all seeds
2691 !
2692       IF (n0.gt.0) THEN
2693       do iters=1,nseed
2694        iseed=is(iters)
2695        index=index+1
2696        movenx(index)=0
2697        parent(1,index)=iseed
2698        parent(2,index)=0
2699
2700        if (vdisulf) then
2701              nss_in(index)=bvar_nss(iseed)
2702              do ij=1,nss_in(index)
2703                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2704                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2705              enddo
2706        endif
2707
2708        do k=1,numch
2709         do j=2,nres-1
2710          do i=1,4
2711           dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2712          enddo
2713         enddo
2714        enddo
2715       enddo
2716       ENDIF
2717 !-----------------------------------------
2718       if (vdisulf) then
2719       do iters=1,nseed
2720        iseed=is(iters)
2721
2722         do k=1,numch
2723           do j=2,nres-1
2724             theta(j+1)=bvar(1,j,k,iseed)
2725             phi(j+2)=bvar(2,j,k,iseed)
2726             alph(j)=bvar(3,j,k,iseed)
2727             omeg(j)=bvar(4,j,k,iseed)
2728           enddo
2729         enddo
2730         call chainbuild
2731
2732 !d       write(iout,*) 'makevar DYNSS',iseed,'#',bvar_ns(iseed),
2733 !d     &     (bvar_s(k,iseed),k=1,bvar_ns(iseed)),
2734 !d     &     bvar_nss(iseed),
2735 !d     &     (bvar_ss(1,k,iseed)-nres,'-',
2736 !d     &      bvar_ss(2,k,iseed)-nres,k=1,bvar_nss(iseed))
2737
2738        do i1=1,bvar_ns(iseed)
2739 !
2740 ! N10 fussion of free halfcysteines in seed 
2741 !      first select CYS with distance < 7A
2742 !
2743          do j1=i1+1,bvar_ns(iseed)
2744            if (dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres) &
2745                .lt.7.0.and. &
2746                iabs(bvar_s(i1,iseed)-bvar_s(j1,iseed)).gt.3) then
2747
2748              index=index+1
2749              movenx(index)=10
2750              parent(1,index)=iseed
2751              parent(2,index)=0
2752              do ij=1,bvar_nss(iseed)
2753                iss_in(ij,index)=bvar_ss(1,ij,iseed)
2754                jss_in(ij,index)=bvar_ss(2,ij,iseed)
2755              enddo
2756              ij=bvar_nss(iseed)+1
2757              nss_in(index)=ij
2758              iss_in(ij,index)=bvar_s(i1,iseed)+nres
2759              jss_in(ij,index)=bvar_s(j1,iseed)+nres
2760
2761 !d             write(iout,*) 'makevar NSS0',index,
2762 !d     &   dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres),
2763 !d     &   nss_in(index),iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres
2764
2765              do k=1,numch
2766               do j=2,nres-1
2767                do i=1,4
2768                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2769                enddo
2770               enddo
2771              enddo
2772
2773            endif
2774          enddo
2775 !
2776 ! N11 type I transdisulfidation
2777 !
2778          do j1=1,bvar_nss(iseed)
2779            if (dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)) &
2780                .lt.7.0.and. &
2781                iabs(bvar_s(i1,iseed)-(bvar_ss(1,j1,iseed)-nres)) &
2782                .gt.3) then
2783
2784              index=index+1
2785              movenx(index)=11
2786              parent(1,index)=iseed
2787              parent(2,index)=0
2788              do ij=1,bvar_nss(iseed)
2789                if (ij.ne.j1) then
2790                 iss_in(ij,index)=bvar_ss(1,ij,iseed)
2791                 jss_in(ij,index)=bvar_ss(2,ij,iseed)
2792                endif
2793              enddo
2794              nss_in(index)=bvar_nss(iseed)
2795              iss_in(j1,index)=bvar_s(i1,iseed)+nres
2796              jss_in(j1,index)=bvar_ss(1,j1,iseed)
2797              if (iss_in(j1,index).gt.jss_in(j1,index)) then
2798                iss_in(j1,index)=bvar_ss(1,j1,iseed)
2799                jss_in(j1,index)=bvar_s(i1,iseed)+nres
2800              endif
2801
2802 !d             write(iout,*) 'makevar NSS1 #1',index,
2803 !d     &       bvar_s(i1,iseed),bvar_ss(1,j1,iseed)-nres,
2804 !d     &       dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)),
2805 !d     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2806 !d     &       ij=1,nss_in(index))
2807
2808              do k=1,numch
2809               do j=2,nres-1
2810                do i=1,4
2811                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2812                enddo
2813               enddo
2814              enddo
2815            endif
2816            if (dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)) &
2817                .lt.7.0.and. &
2818                iabs(bvar_s(i1,iseed)-(bvar_ss(2,j1,iseed)-nres)) &
2819                .gt.3) then
2820
2821              index=index+1
2822              movenx(index)=11
2823              parent(1,index)=iseed
2824              parent(2,index)=0
2825              do ij=1,bvar_nss(iseed)
2826                if (ij.ne.j1) then
2827                 iss_in(ij,index)=bvar_ss(1,ij,iseed)
2828                 jss_in(ij,index)=bvar_ss(2,ij,iseed)
2829                endif
2830              enddo
2831              nss_in(index)=bvar_nss(iseed)
2832              iss_in(j1,index)=bvar_s(i1,iseed)+nres
2833              jss_in(j1,index)=bvar_ss(2,j1,iseed)
2834              if (iss_in(j1,index).gt.jss_in(j1,index)) then
2835                iss_in(j1,index)=bvar_ss(2,j1,iseed)
2836                jss_in(j1,index)=bvar_s(i1,iseed)+nres
2837              endif
2838
2839
2840 !d             write(iout,*) 'makevar NSS1 #2',index,
2841 !d     &       bvar_s(i1,iseed),bvar_ss(2,j1,iseed)-nres,
2842 !d     &       dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)),
2843 !d     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2844 !d     &       ij=1,nss_in(index))
2845
2846              do k=1,numch
2847               do j=2,nres-1
2848                do i=1,4
2849                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2850                enddo
2851               enddo
2852              enddo
2853
2854            endif
2855          enddo
2856        enddo
2857
2858 !
2859 ! N12 type II transdisulfidation
2860 !
2861        do i1=1,bvar_nss(iseed)
2862          do j1=i1+1,bvar_nss(iseed)
2863             if (dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)) &
2864                 .lt.7.0.and. &
2865                 dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)) &
2866                 .lt.7.0.and. &
2867                 iabs(bvar_ss(1,i1,iseed)-bvar_ss(1,j1,iseed)) &
2868                 .gt.3.and. &
2869                 iabs(bvar_ss(2,i1,iseed)-bvar_ss(2,j1,iseed)) &
2870                 .gt.3) then
2871              index=index+1
2872              movenx(index)=12
2873              parent(1,index)=iseed
2874              parent(2,index)=0
2875              do ij=1,bvar_nss(iseed)
2876                if (ij.ne.i1 .and. ij.ne.j1) then
2877                 iss_in(ij,index)=bvar_ss(1,ij,iseed)
2878                 jss_in(ij,index)=bvar_ss(2,ij,iseed)
2879                endif
2880              enddo
2881              nss_in(index)=bvar_nss(iseed)
2882              iss_in(i1,index)=bvar_ss(1,i1,iseed)
2883              jss_in(i1,index)=bvar_ss(1,j1,iseed)
2884              if (iss_in(i1,index).gt.jss_in(i1,index)) then
2885                iss_in(i1,index)=bvar_ss(1,j1,iseed)
2886                jss_in(i1,index)=bvar_ss(1,i1,iseed)
2887              endif
2888              iss_in(j1,index)=bvar_ss(2,i1,iseed)
2889              jss_in(j1,index)=bvar_ss(2,j1,iseed)
2890              if (iss_in(j1,index).gt.jss_in(j1,index)) then
2891                iss_in(j1,index)=bvar_ss(2,j1,iseed)
2892                jss_in(j1,index)=bvar_ss(2,i1,iseed)
2893              endif
2894
2895
2896 !d             write(iout,*) 'makevar NSS2 #1',index,
2897 !d     &       bvar_ss(1,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres, 
2898 !d     &       dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)),
2899 !d     &       bvar_ss(2,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres,
2900 !d     &       dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)),
2901 !d     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2902 !d     &       ij=1,nss_in(index))
2903
2904              do k=1,numch
2905               do j=2,nres-1
2906                do i=1,4
2907                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2908                enddo
2909               enddo
2910              enddo
2911
2912            endif
2913
2914             if (dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)) &
2915                 .lt.7.0.and. &
2916                 dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)) &
2917                 .lt.7.0.and. &
2918                 iabs(bvar_ss(1,i1,iseed)-bvar_ss(2,j1,iseed)) &
2919                 .gt.3.and. &
2920                 iabs(bvar_ss(2,i1,iseed)-bvar_ss(1,j1,iseed)) &
2921                 .gt.3) then
2922              index=index+1
2923              movenx(index)=12
2924              parent(1,index)=iseed
2925              parent(2,index)=0
2926              do ij=1,bvar_nss(iseed)
2927                if (ij.ne.i1 .and. ij.ne.j1) then
2928                 iss_in(ij,index)=bvar_ss(1,ij,iseed)
2929                 jss_in(ij,index)=bvar_ss(2,ij,iseed)
2930                endif
2931              enddo
2932              nss_in(index)=bvar_nss(iseed)
2933              iss_in(i1,index)=bvar_ss(1,i1,iseed)
2934              jss_in(i1,index)=bvar_ss(2,j1,iseed)
2935              if (iss_in(i1,index).gt.jss_in(i1,index)) then
2936                iss_in(i1,index)=bvar_ss(2,j1,iseed)
2937                jss_in(i1,index)=bvar_ss(1,i1,iseed)
2938              endif
2939              iss_in(j1,index)=bvar_ss(2,i1,iseed)
2940              jss_in(j1,index)=bvar_ss(1,j1,iseed)
2941              if (iss_in(j1,index).gt.jss_in(j1,index)) then
2942                iss_in(j1,index)=bvar_ss(1,j1,iseed)
2943                jss_in(j1,index)=bvar_ss(2,i1,iseed)
2944              endif
2945
2946
2947 !d             write(iout,*) 'makevar NSS2 #2',index,
2948 !d     &       bvar_ss(1,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres, 
2949 !d     &       dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)),
2950 !d     &       bvar_ss(2,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres,
2951 !d     &       dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)),
2952 !d     &       (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2953 !d     &       ij=1,nss_in(index))
2954
2955              do k=1,numch
2956               do j=2,nres-1
2957                do i=1,4
2958                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2959                enddo
2960               enddo
2961              enddo
2962
2963            endif
2964
2965
2966          enddo
2967        enddo
2968 !
2969 ! N13 removal of disulfide bond
2970 !
2971       if (bvar_nss(iseed).gt.0) then
2972         i1=bvar_nss(iseed)*ran1(idum)+1
2973
2974              index=index+1
2975              movenx(index)=13
2976              parent(1,index)=iseed
2977              parent(2,index)=0
2978              ij=0
2979              do j1=1,bvar_nss(iseed)
2980               if (j1.ne.i1) then
2981                ij=ij+1
2982                iss_in(ij,index)=bvar_ss(1,j1,iseed)
2983                jss_in(ij,index)=bvar_ss(2,j1,iseed)
2984               endif
2985              enddo
2986              nss_in(index)=bvar_nss(iseed)-1
2987
2988 !d             write(iout,*) 'NSS3',index,i1,
2989 !d     &   bvar_ss(1,i1,iseed)-nres,'=',bvar_ss(2,i1,iseed)-nres,'#',
2990 !d     &   (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres,
2991 !d     &   ij=1,nss_in(index))
2992
2993              do k=1,numch
2994               do j=2,nres-1
2995                do i=1,4
2996                 dihang_in(i,j,k,index)=bvar(i,j,k,iseed)
2997                enddo
2998               enddo
2999              enddo
3000
3001       endif
3002
3003       enddo
3004       endif
3005 !-----------------------------------------
3006
3007
3008
3009       if(index.ne.n) write(iout,*)'make_var : ntry=',index
3010
3011       n=index
3012 !d      do ii=1,n
3013 !d        write (istat,*) "======== ii=",ii," the dihang array"
3014 !d        do i=1,nres
3015 !d       write (istat,'(i5,4f15.5)') i,(dihang_in(k,i,1,ii)*rad2deg,k=1,4)
3016 !d        enddo
3017 !d      enddo
3018       return
3019       end subroutine make_var
3020 !-----------------------------------------------------------------------------
3021       subroutine check_old(icheck,n)
3022
3023 !      implicit real*8 (a-h,o-z)
3024 !      include 'DIMENSIONS'
3025 !      include 'COMMON.CSA'
3026 !      include 'COMMON.BANK'
3027 !      include 'COMMON.CHAIN'
3028 !      include 'COMMON.GEO'
3029       integer :: icheck,n,i1,i2,m,j,i
3030       real(kind=8) :: ctdif,ctdiff,diff,dif
3031
3032       data ctdif /10./
3033       data ctdiff /60./
3034
3035       i1=n
3036       do i2=1,n-1
3037        diff=0.d0
3038        do m=1,numch
3039         do j=2,nres-1
3040          do i=1,4
3041           dif=rad2deg*dabs(dihang_in(i,j,m,i1)-dihang_in(i,j,m,i2))
3042           if(dif.gt.180.0) dif=360.0-dif
3043           if(dif.gt.ctdif) goto 100
3044           diff=diff+dif
3045           if(diff.gt.ctdiff) goto 100
3046          enddo
3047         enddo
3048        enddo
3049        icheck=1
3050        return
3051   100  continue
3052       enddo
3053
3054       icheck=0
3055
3056       return
3057       end subroutine check_old
3058 !-----------------------------------------------------------------------------
3059       subroutine newconf1rr(idum,vvar,nran,i1)
3060
3061 !      implicit real*8 (a-h,o-z)
3062 !      include 'DIMENSIONS'
3063 !      include 'COMMON.IOUNITS'
3064 !      include 'COMMON.CSA'
3065 !      include 'COMMON.BANK'
3066 !      include 'COMMON.CHAIN'
3067 !      include 'COMMON.GEO'
3068 !      real(kind=4) :: ran1,ran2
3069       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3070       integer,dimension(ntotal) :: iold
3071       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3072       real(kind=8) :: ctdif,dif
3073
3074       ctdif=10.
3075
3076       do k=1,numch
3077        do j=2,nres-1
3078         do i=1,4
3079          vvar(i,j,k)=rvar(i,j,k,iseed)
3080         enddo
3081        enddo
3082       enddo
3083
3084       do index=1,nran
3085        iold(index) = 0
3086       enddo
3087
3088       number=ntotgr
3089
3090       iter=0
3091       do index=1,nran
3092    10  iran=  ran1(idum)*number+1
3093        if(iter.gt.number) return
3094        iter=iter+1
3095        if(iter.eq.1) goto 11
3096        do ind=1,index-1
3097         if(iran.eq.iold(ind)) goto 10
3098        enddo
3099    11  continue
3100
3101        do ind=1,ngroup(iran)
3102         i=igroup(1,ind,iran)
3103         j=igroup(2,ind,iran)
3104         k=igroup(3,ind,iran)
3105         dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
3106         if(dif.gt.180.) dif=360.-dif
3107         if(dif.gt.ctdif) goto 20
3108        enddo
3109        if(iter.gt.number) goto 20
3110        goto 10
3111    20  continue
3112        do ind=1,ngroup(iran)
3113         i=igroup(1,ind,iran)
3114         j=igroup(2,ind,iran)
3115         k=igroup(3,ind,iran)
3116         vvar(i,j,k)=rvar(i,j,k,i1)
3117        enddo
3118        iold(index)=iran
3119       enddo
3120
3121       return
3122       end subroutine newconf1rr
3123 !-----------------------------------------------------------------------------
3124       subroutine newconf1br(idum,vvar,nran,i1)
3125
3126       use energy_data, only: ndih_nconstr,idih_nconstr
3127       use control_data, only: i2ndstr
3128 !      implicit real*8 (a-h,o-z)
3129 !      include 'DIMENSIONS'
3130 !      include 'COMMON.IOUNITS'
3131 !      include 'COMMON.CSA'
3132 !      include 'COMMON.BANK'
3133 !      include 'COMMON.CHAIN'
3134 !      include 'COMMON.GEO'
3135 !      include 'COMMON.TORCNSTR'
3136 !      include 'COMMON.CONTROL'
3137 !      real(kind=4) :: ran1,ran2
3138       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3139       integer,dimension(ntotal) :: iold
3140       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,juhc,ind
3141       real(kind=8) :: ctdif,dif,rtmp
3142
3143       ctdif=10.
3144
3145       do k=1,numch
3146        do j=2,nres-1
3147         do i=1,4
3148          vvar(i,j,k)=bvar(i,j,k,iseed)
3149         enddo
3150        enddo
3151       enddo
3152
3153       do index=1,nran
3154        iold(index) = 0
3155       enddo
3156
3157       number=ntotgr
3158
3159       iter=0
3160       do index=1,nran
3161    10  iran=  ran1(idum)*number+1
3162        if(i2ndstr.gt.0) then
3163         rtmp=ran1(idum)
3164         if(rtmp.le.rdih_bias) then
3165          i=0
3166          do j=1,ndih_nconstr
3167           if(igroup(2,1,iran).eq.idih_nconstr(j))i=j
3168          enddo
3169          if(i.eq.0) then
3170           juhc=0
3171 4321      juhc=juhc+1
3172           iran=  ran1(idum)*number+1
3173           i=0
3174           do j=1,ndih_nconstr
3175            if(igroup(2,1,iran).eq.idih_nconstr(j))i=j
3176           enddo
3177           if(i.eq.0.or.juhc.lt.1000)goto 4321
3178           if(juhc.eq.1000) then
3179            print *, 'move 6 : failed to find unconstrained group'
3180            write(iout,*) 'move 6 : failed to find unconstrained group'
3181           endif
3182          endif
3183         endif
3184        endif
3185        if(iter.gt.number) return
3186        iter=iter+1
3187        if(iter.eq.1) goto 11
3188        do ind=1,index-1
3189         if(iran.eq.iold(ind)) goto 10
3190        enddo
3191    11  continue
3192
3193        do ind=1,ngroup(iran)
3194         i=igroup(1,ind,iran)
3195         j=igroup(2,ind,iran)
3196         k=igroup(3,ind,iran)
3197         dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
3198         if(dif.gt.180.) dif=360.-dif
3199         if(dif.gt.ctdif) goto 20
3200        enddo
3201        if(iter.gt.number) goto 20
3202        goto 10
3203    20  continue
3204        do ind=1,ngroup(iran)
3205         i=igroup(1,ind,iran)
3206         j=igroup(2,ind,iran)
3207         k=igroup(3,ind,iran)
3208         vvar(i,j,k)=rvar(i,j,k,i1)
3209        enddo
3210        iold(index)=iran
3211       enddo
3212
3213       return
3214       end subroutine newconf1br
3215 !-----------------------------------------------------------------------------
3216       subroutine newconf1bb(idum,vvar,nran,i1)
3217
3218 !      implicit real*8 (a-h,o-z)
3219 !      include 'DIMENSIONS'
3220 !      include 'COMMON.IOUNITS'
3221 !      include 'COMMON.CSA'
3222 !      include 'COMMON.BANK'
3223 !      include 'COMMON.CHAIN'
3224 !      include 'COMMON.GEO'
3225 !      real(kind=4) :: ran1,ran2
3226       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3227       integer,dimension(ntotal) :: iold
3228       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3229       real(kind=8) :: ctdif,dif
3230
3231       ctdif=10.
3232
3233       do k=1,numch
3234        do j=2,nres-1
3235         do i=1,4
3236          vvar(i,j,k)=bvar(i,j,k,iseed)
3237         enddo
3238        enddo
3239       enddo
3240
3241       do index=1,nran
3242        iold(index) = 0
3243       enddo
3244
3245       number=ntotgr
3246
3247       iter=0
3248       do index=1,nran
3249    10  iran=  ran1(idum)*number+1
3250        if(iter.gt.number) return
3251        iter=iter+1
3252        if(iter.eq.1) goto 11
3253        do ind=1,index-1
3254         if(iran.eq.iold(ind)) goto 10
3255        enddo
3256    11  continue
3257
3258        do ind=1,ngroup(iran)
3259         i=igroup(1,ind,iran)
3260         j=igroup(2,ind,iran)
3261         k=igroup(3,ind,iran)
3262         dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
3263         if(dif.gt.180.) dif=360.-dif
3264         if(dif.gt.ctdif) goto 20
3265        enddo
3266        if(iter.gt.number) goto 20
3267        goto 10
3268    20  continue
3269        do ind=1,ngroup(iran)
3270         i=igroup(1,ind,iran)
3271         j=igroup(2,ind,iran)
3272         k=igroup(3,ind,iran)
3273         vvar(i,j,k)=bvar(i,j,k,i1)
3274        enddo
3275        iold(index)=iran
3276       enddo
3277
3278       return
3279       end subroutine newconf1bb
3280 !-----------------------------------------------------------------------------
3281       subroutine newconf1arr(idum,vvar,nran,i1)
3282
3283 !      implicit real*8 (a-h,o-z)
3284 !      include 'DIMENSIONS'
3285 !      include 'COMMON.IOUNITS'
3286 !      include 'COMMON.CSA'
3287 !      include 'COMMON.BANK'
3288 !      include 'COMMON.CHAIN'
3289 !      include 'COMMON.GEO'
3290 !      real(kind=4) :: ran1,ran2
3291       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3292       integer,dimension(ntotal) :: iold
3293       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3294       real(kind=8) :: ctdif,dif
3295
3296       ctdif=10.
3297
3298       do k=1,numch
3299        do j=2,nres-1
3300         do i=1,4
3301          vvar(i,j,k)=rvar(i,j,k,iseed)
3302         enddo
3303        enddo
3304       enddo
3305
3306       do index=1,nran
3307        iold(index) = 0
3308       enddo
3309
3310       number=nres-2
3311
3312       iter=0
3313       do index=1,nran
3314    10  iran=  ran1(idum)*number+1
3315        if(iter.gt.number) return
3316        iter=iter+1
3317        if(iter.eq.1) goto 11
3318        do ind=1,index-1
3319         if(iran.eq.iold(ind)) goto 10
3320        enddo
3321    11  continue
3322
3323        do ind=1,ngroup(iran)
3324         i=igroup(1,ind,iran)
3325         j=igroup(2,ind,iran)
3326         k=igroup(3,ind,iran)
3327         dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
3328         if(dif.gt.180.) dif=360.-dif
3329         if(dif.gt.ctdif) goto 20
3330        enddo
3331        if(iter.gt.number) goto 20
3332        goto 10
3333    20  continue
3334        do ind=1,ngroup(iran)
3335         i=igroup(1,ind,iran)
3336         j=igroup(2,ind,iran)
3337         k=igroup(3,ind,iran)
3338         vvar(i,j,k)=rvar(i,j,k,i1)
3339        enddo
3340        iold(index)=iran
3341       enddo
3342
3343       return
3344       end subroutine newconf1arr
3345 !-----------------------------------------------------------------------------
3346       subroutine newconf1abr(idum,vvar,nran,i1)
3347
3348       use energy_data, only: ndih_nconstr,idih_nconstr
3349       use control_data, only: i2ndstr
3350 !      implicit real*8 (a-h,o-z)
3351 !      include 'DIMENSIONS'
3352 !      include 'COMMON.IOUNITS'
3353 !      include 'COMMON.CSA'
3354 !      include 'COMMON.BANK'
3355 !      include 'COMMON.CHAIN'
3356 !      include 'COMMON.GEO'
3357 !      include 'COMMON.TORCNSTR'
3358 !      include 'COMMON.CONTROL'
3359 !      real(kind=4) :: ran1,ran2
3360       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3361       integer,dimension(ntotal) :: iold
3362       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3363       real(kind=8) :: ctdif,dif,rtmp
3364
3365       ctdif=10.
3366
3367       do k=1,numch
3368        do j=2,nres-1
3369         do i=1,4
3370          vvar(i,j,k)=bvar(i,j,k,iseed)
3371         enddo
3372        enddo
3373       enddo
3374
3375       do index=1,nran
3376        iold(index) = 0
3377       enddo
3378
3379       number=nres-2
3380
3381       iter=0
3382       do index=1,nran
3383    10  iran=  ran1(idum)*number+1
3384       if(i2ndstr.gt.0) then
3385        rtmp=ran1(idum)
3386        if(rtmp.le.rdih_bias) then
3387         iran=ran1(idum)*ndih_nconstr+1
3388         iran=idih_nconstr(iran)
3389        endif
3390       endif
3391        if(iter.gt.number) return
3392        iter=iter+1
3393        if(iter.eq.1) goto 11
3394        do ind=1,index-1
3395         if(iran.eq.iold(ind)) goto 10
3396        enddo
3397    11  continue
3398
3399        do ind=1,ngroup(iran)
3400         i=igroup(1,ind,iran)
3401         j=igroup(2,ind,iran)
3402         k=igroup(3,ind,iran)
3403         dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1))
3404         if(dif.gt.180.) dif=360.-dif
3405         if(dif.gt.ctdif) goto 20
3406        enddo
3407        if(iter.gt.number) goto 20
3408        goto 10
3409    20  continue
3410        do ind=1,ngroup(iran)
3411         i=igroup(1,ind,iran)
3412         j=igroup(2,ind,iran)
3413         k=igroup(3,ind,iran)
3414         vvar(i,j,k)=rvar(i,j,k,i1)
3415        enddo
3416        iold(index)=iran
3417       enddo
3418
3419       return
3420       end subroutine newconf1abr
3421 !-----------------------------------------------------------------------------
3422       subroutine newconf1abb(idum,vvar,nran,i1)
3423
3424       use energy_data, only: ndih_nconstr,idih_nconstr
3425       use control_data, only: i2ndstr
3426 !      implicit real*8 (a-h,o-z)
3427 !      include 'DIMENSIONS'
3428 !      include 'COMMON.IOUNITS'
3429 !      include 'COMMON.CSA'
3430 !      include 'COMMON.BANK'
3431 !      include 'COMMON.CHAIN'
3432 !      include 'COMMON.GEO'
3433 !      include 'COMMON.TORCNSTR'
3434 !      include 'COMMON.CONTROL'
3435 !      real(kind=4) :: ran1,ran2
3436       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3437       integer,dimension(ntotal) :: iold
3438       integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind
3439       real(kind=8) :: ctdif,dif,rtmp
3440
3441       ctdif=10.
3442
3443       do k=1,numch
3444        do j=2,nres-1
3445         do i=1,4
3446          vvar(i,j,k)=bvar(i,j,k,iseed)
3447         enddo
3448        enddo
3449       enddo
3450
3451       do index=1,nran
3452        iold(index) = 0
3453       enddo
3454
3455       number=nres-2
3456
3457       iter=0
3458       do index=1,nran
3459    10  iran=  ran1(idum)*number+1
3460       if(i2ndstr.gt.0) then
3461        rtmp=ran1(idum)
3462        if(rtmp.le.rdih_bias) then
3463         iran=ran1(idum)*ndih_nconstr+1
3464         iran=idih_nconstr(iran)
3465        endif
3466       endif
3467        if(iter.gt.number) return
3468        iter=iter+1
3469        if(iter.eq.1) goto 11
3470        do ind=1,index-1
3471         if(iran.eq.iold(ind)) goto 10
3472        enddo
3473    11  continue
3474
3475        do ind=1,ngroup(iran)
3476         i=igroup(1,ind,iran)
3477         j=igroup(2,ind,iran)
3478         k=igroup(3,ind,iran)
3479         dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
3480         if(dif.gt.180.) dif=360.-dif
3481         if(dif.gt.ctdif) goto 20
3482        enddo
3483        if(iter.gt.number) goto 20
3484        goto 10
3485    20  continue
3486        do ind=1,ngroup(iran)
3487         i=igroup(1,ind,iran)
3488         j=igroup(2,ind,iran)
3489         k=igroup(3,ind,iran)
3490         vvar(i,j,k)=bvar(i,j,k,i1)
3491        enddo
3492        iold(index)=iran
3493       enddo
3494
3495       return
3496       end subroutine newconf1abb
3497 !-----------------------------------------------------------------------------
3498       subroutine newconf_residue(idum,vvar,i1,isize)
3499
3500       use energy_data, only: ndih_nconstr,idih_nconstr
3501       use control_data, only: i2ndstr
3502       use MPI_data
3503       include 'mpif.h'
3504 !      implicit real*8 (a-h,o-z)
3505 !      include 'DIMENSIONS'
3506 !      include 'COMMON.IOUNITS'
3507 !      include 'COMMON.CSA'
3508 !      include 'COMMON.BANK'
3509 !      include 'COMMON.CHAIN'
3510 !      include 'COMMON.GEO'
3511 !      include 'COMMON.TORCNSTR'
3512 !      include 'COMMON.CONTROL'
3513 !      real(kind=4) :: ran1,ran2
3514       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3515       integer,dimension(ntotal) :: iold
3516       integer :: i,j,k,idum,i1,isize,iran,number,iter,ind,iend,istart,&
3517             ierror,ierrcode
3518       real(kind=8) :: ctdif,dif,rtmp
3519
3520       ctdif=10.
3521
3522       if (iseed.gt.mxio .or. iseed.lt.1) then
3523         write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
3524         call mpi_abort(mpi_comm_world,ierror,ierrcode)
3525       endif
3526       do k=1,numch
3527        do j=2,nres-1
3528         do i=1,4
3529          vvar(i,j,k)=bvar(i,j,k,iseed)
3530         enddo
3531        enddo
3532       enddo
3533
3534
3535       k=1
3536       number=nres+isize-2
3537       iter=1
3538    10 iran=  ran1(idum)*number+1
3539       if(i2ndstr.gt.0) then
3540        rtmp=ran1(idum)
3541        if(rtmp.le.rdih_bias) then
3542         iran=ran1(idum)*ndih_nconstr+1
3543         iran=idih_nconstr(iran)
3544        endif
3545       endif
3546       istart=iran-isize+1
3547       iend=iran
3548       if(istart.lt.2) istart=2
3549       if(iend.gt.nres-1) iend=nres-1
3550
3551       if(iter.eq.1) goto 11
3552       do ind=1,iter-1
3553        if(iran.eq.iold(ind)) goto 10
3554       enddo
3555    11 continue
3556
3557       do j=istart,iend
3558        do i=1,4
3559          dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
3560          if(dif.gt.180.) dif=360.-dif
3561          if(dif.gt.ctdif) goto 20
3562        enddo
3563       enddo
3564       iold(iter)=iran
3565       iter=iter+1
3566       if(iter.gt.number) goto 20
3567       goto 10
3568
3569    20 continue
3570       do j=istart,iend
3571        do i=1,4
3572         vvar(i,j,k)=bvar(i,j,k,i1)
3573        enddo
3574       enddo
3575
3576       return
3577       end subroutine newconf_residue
3578 !-----------------------------------------------------------------------------
3579       subroutine newconf_copy(idum,vvar,i1,istart,iend)
3580
3581       use MPI_data
3582       include 'mpif.h'
3583 !      implicit real*8 (a-h,o-z)
3584 !      include 'DIMENSIONS'
3585 !      include 'COMMON.IOUNITS'
3586 !      include 'COMMON.CSA'
3587 !      include 'COMMON.BANK'
3588 !      include 'COMMON.CHAIN'
3589 !      include 'COMMON.GEO'
3590 !      include 'COMMON.TORCNSTR'
3591 !      include 'COMMON.CONTROL'
3592 !      real(kind=4) :: ran1,ran2
3593       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3594       integer,dimension(ntotal) :: iold
3595       integer :: i,j,k,idum,i1,istart,iend,ierror,ierrcode
3596       real(kind=8) :: ctdif,dif
3597
3598       ctdif=10.
3599
3600       if (iseed.gt.mxio .or. iseed.lt.1) then
3601         write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
3602         call mpi_abort(mpi_comm_world,ierror,ierrcode)
3603       endif
3604       do k=1,numch
3605        do j=2,nres-1
3606         do i=1,4
3607          vvar(i,j,k)=bvar(i,j,k,iseed)
3608         enddo
3609        enddo
3610       enddo
3611
3612
3613       do j=istart,iend
3614        do i=1,4
3615         vvar(i,j,1)=bvar(i,j,1,i1)
3616        enddo
3617       enddo
3618
3619       return
3620       end subroutine newconf_copy
3621 !-----------------------------------------------------------------------------
3622       subroutine newconf_residue_hairpin(idum,vvar,i1,fail)
3623
3624       use geometry_data
3625 !      use random, only: iran_num
3626       use MPI_data
3627       use compare, only:hairpin
3628      
3629       include 'mpif.h'
3630 !      implicit real*8 (a-h,o-z)
3631 !      include 'DIMENSIONS'
3632 !      include 'COMMON.IOUNITS'
3633 !      include 'COMMON.CSA'
3634 !      include 'COMMON.BANK'
3635 !      include 'COMMON.CHAIN'
3636 !      include 'COMMON.GEO'
3637 !      include 'COMMON.VAR'
3638 !      real(kind=4) :: ran1,ran2
3639       real(kind=8),dimension(mxang,nres,mxch) :: vvar   !(mxang,maxres,mxch)
3640       integer,dimension(ntotal) :: iold
3641       integer :: nharp,iharp(4,nres/3),icipa(nres/3)
3642       logical :: fail,not_done
3643       integer :: idum,i,j,k,i1,iend,istart,iii,n_used,icount,iih,&
3644             ierror,ierrcode
3645       real(kind=8) :: ctdif,dif
3646
3647       ctdif=10.
3648
3649       fail=.false.
3650       if (iseed.gt.mxio .or. iseed.lt.1) then
3651         write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed
3652         call mpi_abort(mpi_comm_world,ierror,ierrcode)
3653       endif
3654       do k=1,numch
3655        do j=2,nres-1
3656         do i=1,4
3657          vvar(i,j,k)=bvar(i,j,k,iseed)
3658         enddo
3659        enddo
3660       enddo
3661       do k=1,numch
3662       do j=2,nres-1
3663         theta(j+1)=bvar(1,j,k,i1)
3664         phi(j+2)=bvar(2,j,k,i1)
3665         alph(j)=bvar(3,j,k,i1)
3666         omeg(j)=bvar(4,j,k,i1)
3667       enddo
3668       enddo
3669 !      call intout
3670       call chainbuild
3671       call hairpin(.false.,nharp,iharp)
3672
3673       if (nharp.eq.0) then
3674         fail=.true.
3675         return
3676       endif
3677
3678       n_used=0
3679
3680       DO III=1,NHARP
3681
3682       not_done = .true.
3683       icount=0
3684       do while (not_done)
3685         icount=icount+1
3686         iih=iran_num(1,nharp)
3687         do k=1,n_used
3688           if (iih.eq.icipa(k)) then
3689             iih=0
3690             goto 22
3691           endif
3692         enddo
3693         not_done=.false.
3694         n_used=n_used+1
3695         icipa(n_used)=iih
3696    22   continue
3697         not_done = not_done .and. icount.le.nharp
3698       enddo
3699
3700       if (iih.eq.0) then
3701         write (iout,*) "CHUJ NASTAPIL W NEWCONF_RESIDUE_HAIRPIN!!!!"
3702         fail=.true.
3703         return
3704       endif
3705       
3706       istart=iharp(1,iih)+1
3707       iend=iharp(2,iih)
3708
3709 !dd      write (iout,*) "newconf_residue_hairpin: iih",iih,
3710 !dd     &   " istart",istart," iend",iend
3711
3712       do k=1,numch
3713       do j=istart,iend
3714        do i=1,4
3715          dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1))
3716          if(dif.gt.180.) dif=360.-dif
3717          if(dif.gt.ctdif) goto 20
3718        enddo
3719       enddo
3720       enddo
3721       goto 10
3722    20 continue
3723       do k=1,numch
3724       do j=istart,iend
3725        do i=1,4
3726         vvar(i,j,k)=bvar(i,j,k,i1)
3727        enddo
3728       enddo
3729       enddo
3730 !      do j=1,numch
3731 !       do l=2,nres-1
3732 !        write (iout,'(4f8.3)') (rad2deg*vvar(i,l,j),i=1,4)
3733 !       enddo
3734 !      enddo
3735       return
3736    10 continue
3737       ENDDO
3738
3739       fail=.true.
3740
3741       return
3742       end subroutine newconf_residue_hairpin
3743 !-----------------------------------------------------------------------------
3744       subroutine gen_hairpin
3745
3746       use geometry_data
3747       use MD_data
3748       use compare, only:hairpin
3749 !      implicit real*8 (a-h,o-z)
3750 !      include 'DIMENSIONS'
3751 !      include 'COMMON.IOUNITS'
3752 !      include 'COMMON.CSA'
3753 !      include 'COMMON.BANK'
3754 !      include 'COMMON.CHAIN'
3755 !      include 'COMMON.GEO'
3756 !      include 'COMMON.VAR'
3757 !      include 'COMMON.HAIRPIN'
3758       integer :: i1,j,k,iters
3759
3760 !      write (iout,*) 'Entering GEN_HAIRPIN'
3761       do iters=1,nseed
3762         i1=is(iters)
3763         do k=1,numch
3764           do j=2,nres-1
3765             theta(j+1)=bvar(1,j,k,i1)
3766             phi(j+2)=bvar(2,j,k,i1)
3767             alph(j)=bvar(3,j,k,i1)
3768             omeg(j)=bvar(4,j,k,i1)
3769           enddo
3770         enddo
3771         call chainbuild
3772         call hairpin(.false.,nharp_seed(iters),iharp_seed(1,1,iters))
3773       enddo
3774
3775       nharp_tot=0
3776       do iters=1,nseed
3777         nharp_tot=nharp_tot+nharp_seed(iters)
3778         nharp_use(iters)=4*nharp_seed(iters)
3779         do j=1,nharp_seed(iters)
3780           iharp_use(0,j,iters)=4
3781           do k=1,4
3782             iharp_use(k,j,iters)=0
3783           enddo
3784         enddo
3785       enddo
3786
3787       write (iout,*) 'GEN_HAIRPIN: nharp_tot',nharp_tot
3788 !dd      do i=1,nseed
3789 !dd        write (iout,*) 'seed',i 
3790 !dd        write (iout,*) 'nharp_seed',nharp_seed(i),
3791 !dd     &     ' nharp_use',nharp_use(i)
3792 !d        write (iout,*) 'iharp_seed, iharp_use'
3793 !d        do j=1,nharp_seed(i)
3794 !d          write (iout,'(7i3)') iharp_seed(1,j,i),iharp_seed(2,j,i),
3795 !d     &      (iharp_use(k,j,i),k=0,4) 
3796 !d        enddo
3797 !dd      enddo
3798       return
3799       end subroutine gen_hairpin
3800 !-----------------------------------------------------------------------------
3801       subroutine select_frag(nn,nh,nl,ns,nb,i_csa)
3802
3803       use geometry_data
3804       use MD_data
3805       use compare_data
3806 !      implicit real*8 (a-h,o-z)
3807 !      include 'DIMENSIONS'
3808 !      include 'COMMON.IOUNITS'
3809 !      include 'COMMON.CSA'
3810 !      include 'COMMON.BANK'
3811 !      include 'COMMON.CHAIN'
3812 !      include 'COMMON.GEO'
3813 !      include 'COMMON.VAR'
3814 !      include 'COMMON.HAIRPIN'
3815 !      include 'COMMON.DISTFIT'
3816       character(len=50) :: linia
3817       integer :: isec(nres)
3818       integer :: i,j,i1,k,nn,nh,nl,ns,nb,i_csa,nl1,ns1
3819
3820       nn=0
3821       nh=0
3822       nl=0
3823       ns=0
3824       nb=0
3825 !d      write (iout,*) 'Entering select_frag'
3826       do i1=1,nbank
3827         do i=1,nres
3828           isec(i)=0
3829         enddo
3830         do k=1,numch
3831           do j=2,nres-1
3832             theta(j+1)=bvar(1,j,k,i1)
3833             phi(j+2)=bvar(2,j,k,i1)
3834             alph(j)=bvar(3,j,k,i1)
3835             omeg(j)=bvar(4,j,k,i1)
3836           enddo
3837         enddo
3838         call chainbuild
3839 !d        write (iout,*) ' -- ',i1,' -- ' 
3840         call secondary2(.false.)
3841 !
3842 ! bvar_frag nn==pair of nonlocal strands in beta sheet (loop>4)
3843 !               strands > 4 residues; used by N7 and N16
3844 !
3845         do j=1,nbfrag
3846 !
3847 !test 09/12/02 bfrag(2,j)-bfrag(1,j).gt.3
3848 !
3849               do i=bfrag(1,j),bfrag(2,j)
3850                 isec(i)=1
3851               enddo
3852               do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
3853                 isec(i)=1
3854               enddo
3855
3856            if ( (bfrag(3,j).lt.bfrag(4,j) .or. &
3857                 bfrag(4,j)-bfrag(2,j).gt.4) .and. &
3858                 bfrag(2,j)-bfrag(1,j).gt.4 ) then
3859              nn=nn+1
3860
3861
3862              if (bfrag(3,j).lt.bfrag(4,j)) then
3863                write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
3864                            "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
3865                                 ",",bfrag(3,j)-1,"-",bfrag(4,j)-1
3866              else
3867                write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') &
3868                            "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,&
3869                                 ",",bfrag(4,j)-1,"-",bfrag(3,j)-1
3870
3871              endif
3872 !d             call write_pdb(i_csa*1000+nn+nh,linia,0d0)
3873
3874              bvar_frag(nn,1)=i1
3875              bvar_frag(nn,2)=4
3876              do i=1,4
3877                bvar_frag(nn,i+2)=bfrag(i,j)
3878              enddo
3879            endif
3880         enddo
3881
3882 !
3883 ! hvar_frag nh==helices; used by N8 and N9
3884 !
3885         do j=1,nhfrag
3886
3887               do i=hfrag(1,j),hfrag(2,j)
3888                 isec(i)=2
3889               enddo
3890
3891           if ( hfrag(2,j)-hfrag(1,j).gt.4 ) then
3892             nh=nh+1
3893
3894 !d            write(linia,'(a6,i3,a1,i3)')
3895 !d     &                "select",hfrag(1,j)-1,"-",hfrag(2,j)-1
3896 !d            call write_pdb(i_csa*1000+nn+nh,linia,0d0)
3897
3898             hvar_frag(nh,1)=i1
3899             hvar_frag(nh,2)=hfrag(1,j)
3900             hvar_frag(nh,3)=hfrag(2,j)
3901           endif
3902         enddo
3903
3904         
3905 !v        write(iout,'(i4,1pe12.4,1x,1000i1)') 
3906 !v     &         i1,bene(i1),(isec(i),i=1,nres)
3907 !v            write(linia,'(i4,1x,1000i1)') 
3908 !v     &         i1,(isec(i),i=1,nres)
3909 !v            call write_pdb(i_csa*1000+i1,linia,bene(i1))
3910 !
3911 ! lvar_frag nl==loops; used by N14
3912 !
3913         i=1
3914         nl1=nl
3915         do while (i.lt.nres)
3916           if (isec(i).eq.0) then
3917            nl=nl+1
3918            lvar_frag(nl,1)=i1
3919            lvar_frag(nl,2)=i
3920            i=i+1
3921            do while (isec(i).eq.0.and.i.le.nres)
3922              i=i+1
3923            enddo
3924            lvar_frag(nl,3)=i-1
3925            if (lvar_frag(nl,3)-lvar_frag(nl,2).lt.1) nl=nl-1
3926           endif
3927           i=i+1
3928         enddo
3929 !d        write(iout,'(4i5)') (i,(lvar_frag(i,ii),ii=1,3),i=nl1+1,nl)
3930
3931 !
3932 ! svar_frag ns==an secondary structure element; used by N15
3933 !
3934         i=1
3935         ns1=ns
3936         do while (i.lt.nres)
3937           if (isec(i).gt.0) then
3938            ns=ns+1
3939            svar_frag(ns,1)=i1
3940            svar_frag(ns,2)=i
3941            i=i+1
3942            do while (isec(i).gt.0.and.isec(i-1).eq.isec(i) &
3943                                .and.i.le.nres)
3944              i=i+1
3945            enddo
3946            svar_frag(ns,3)=i-1
3947            if (svar_frag(ns,3)-svar_frag(ns,2).lt.1) ns=ns-1
3948           endif
3949           if (isec(i).eq.0) i=i+1
3950         enddo
3951 !d        write(iout,'(4i5)') (i,(svar_frag(i,ii),ii=1,3),i=ns1+1,ns)
3952
3953 !
3954 ! avar_frag nb==any pair of beta strands; used by N17
3955 !
3956         do j=1,nbfrag
3957              nb=nb+1
3958              avar_frag(nb,1)=i1
3959              do i=1,4
3960                avar_frag(nb,i+1)=bfrag(i,j)
3961              enddo
3962         enddo
3963
3964       enddo
3965
3966       return
3967       end subroutine select_frag
3968 !-----------------------------------------------------------------------------
3969 ! together.F
3970 !-----------------------------------------------------------------------------
3971       subroutine together
3972
3973 !  feeds tasks for parallel processing
3974       use MPI_data
3975       use geometry_data
3976       use control_data, only: vdisulf
3977       use energy_data
3978       use io, only:from_int,write_csa_pdb
3979 !      implicit real*8 (a-h,o-z)
3980 !      include 'DIMENSIONS'
3981       include 'mpif.h'
3982 !      real(kind=4) :: ran1,ran2
3983 !      include 'COMMON.CSA'
3984 !      include 'COMMON.BANK'
3985 !      include 'COMMON.IOUNITS'
3986 !      include 'COMMON.CHAIN'
3987 !      include 'COMMON.TIME1'
3988 !      include 'COMMON.SETUP'
3989 !      include 'COMMON.VAR'
3990 !      include 'COMMON.GEO'
3991 !      include 'COMMON.CONTROL'
3992 !      include 'COMMON.SBRIDGE'
3993       real(kind=4) :: tcpu
3994       real(kind=8) :: time_start,time_start_c,time0f,time0i
3995       logical :: ovrtim,sync_iter,timeout,flag,timeout1
3996       integer,dimension(mpi_status_size) :: muster
3997       real(kind=8),dimension(0:100) :: t100
3998       integer,dimension(mxio) :: indx
3999       real(kind=8),dimension(6*nres) :: xout    !(maxvar) (maxvar=6*maxres)
4000       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
4001       integer,dimension(9) :: ind
4002       real(kind=8),dimension(2) :: cout
4003       real(kind=8),parameter :: rad=1.745329252d-2
4004
4005       integer :: i,m,j,jlee,nft,idum,nrmsdb,nrmsdb1,ierr,ierror,ierrcode,&
4006             ntrial,ntry,idum2,imax,idumm,nconfr,iconf,mm,k,im,nst,ifar,&
4007             iter,irecv,isent,iw_pdb,nft0i,nft00_c,nft00,ifrom,ij,&
4008             ireq,ireq2
4009       real(kind=8) :: adif,p_cut,cutdifr,rmsdbc1c,time1i,ctdif1,xctdif,&
4010                  time2i,tstart,tend1
4011 !ccccccccccccccccccccccccccccccccccccccccccccccc
4012       sync_iter=.true.   !el
4013       nft=0     !el
4014       time_start=0.0d0
4015       IF (ME.EQ.KING) THEN
4016
4017        time0f=MPI_WTIME()
4018        ilastnstep=1
4019        sync_iter=.false.
4020        numch=1
4021        nrmsdb=0
4022        nrmsdb1=0
4023        rmsdbc1c=rmsdbc1
4024        nstep=0
4025        call csa_read
4026        call make_array
4027
4028        if(iref.ne.0) call from_int(1,0,idum)
4029
4030 ! To minimize input conformation (bank conformation)
4031 ! Output to $mol.reminimized
4032        if (irestart.lt.0) then
4033         call read_bank(0,nft,cutdifr)
4034         if (irestart.lt.-10) then
4035          p_cut=nres*4.d0
4036          call prune_bank(p_cut)
4037          return
4038         endif
4039         call reminimize(jlee)
4040         return
4041        endif
4042
4043        if (irestart.eq.0) then
4044         call initial_write
4045         nbank=nconf
4046         ntbank=nconf
4047         if (ntbankm.eq.0) ntbank=0
4048         nstep=0
4049         nft=0
4050         do i=1,mxio
4051          ibank(i)=0
4052          jbank(i)=0
4053         enddo
4054        else
4055         call restart_write
4056 !!bankt        call read_bankt(jlee,nft,cutdifr)
4057         call read_bank(jlee,nft,cutdifr)
4058         call read_rbank(jlee,adif)
4059         if(iref.ne.0) call from_int(1,0,idum)
4060        endif
4061
4062        nstmax=nstmax+nstep
4063        ntrial=n1+n2+n3+n4+n5+n6+n7+n8
4064        ntry=ntrial+1
4065        ntry=ntry*nseed
4066
4067 ! ntrial : number of trial conformations per seed.
4068 ! ntry : total number of trial conformations including seed conformations.
4069
4070        idum2=-123
4071 !       imax=2**31-1
4072        imax=huge(0)
4073        ENDIF
4074
4075        call mpi_bcast(jend,1,mpi_integer,0,CG_COMM,ierr)
4076 !ccccccccccccccccccccccccccccccccccccccc
4077        do 300 jlee=1,jend
4078 !ccccccccccccccccccccccccccccccccccccccc
4079   331   continue 
4080         IF (ME.EQ.KING) THEN
4081         if(sync_iter) goto 333
4082         idum=-  ran2(idum2)*imax
4083         if(jlee.lt.jstart) goto 300
4084
4085 ! Restart the random number generator for conformation generation
4086
4087         if(irestart.gt.0) then
4088          idum2=idum2+nstep
4089          if(idum2.le.0) idum2=-idum2+1
4090          idum=-  ran2(idum2)*imax
4091         endif
4092
4093         idumm=idum
4094         call vrndst(idumm)
4095
4096         open(icsa_seed,file=csa_seed,status="old")
4097         write(icsa_seed,*) "jlee : ",jlee
4098         close(icsa_seed)
4099
4100       call history_append
4101       write(icsa_history,*) "number of procs is ",nodes
4102       write(icsa_history,*) jlee,idum,idum2
4103       close(icsa_history)
4104
4105 !ccccccccccccccccccccccccccccccccccccccccccccccc
4106   333 icycle=0
4107
4108        call history_append
4109         write(icsa_history,*) "nbank is ",nbank
4110        close(icsa_history)
4111
4112       if(irestart.eq.1) goto 111
4113       if(irestart.eq.2) then
4114        icycle=0
4115        do i=1,nbank
4116         ibank(i)=1
4117        enddo
4118        do i=nbank+1,nbank+nconf
4119         ibank(i)=0
4120        enddo
4121       endif
4122
4123 !  start energy minimization
4124       nconfr=max0(nconf+nadd,nodes-1)
4125       if (sync_iter) nconf_in=0
4126 !  king-emperor - feed input and sort output
4127        write (iout,*) "NCONF_IN",nconf_in
4128        m=0
4129        if (nconf_in.gt.0) then
4130 ! al 7/2/00 - added possibility to read in some of the initial conformations
4131         do m=1,nconf_in
4132           read (intin,'(i5)',end=11,err=12) iconf
4133    12     continue
4134           write (iout,*) "write READ_ANGLES",iconf,m
4135           call read_angles(intin,*11)
4136           if (iref.eq.0) then
4137             mm=m
4138           else
4139             mm=m+1
4140           endif
4141           do j=2,nres-1
4142             dihang_in(1,j,1,mm)=theta(j+1)
4143             dihang_in(2,j,1,mm)=phi(j+2)
4144             dihang_in(3,j,1,mm)=alph(j)
4145             dihang_in(4,j,1,mm)=omeg(j)
4146           enddo
4147         enddo ! m
4148         goto 13
4149    11   write (iout,*) nconf_in," conformations requested, but only",&
4150          m-1," found in the angle file."
4151         nconf_in=m-1
4152    13   continue
4153         m=nconf_in
4154         write (iout,*) nconf_in,&
4155           " initial conformations have been read in."
4156        endif
4157        if (iref.eq.0) then
4158         if (nconfr.gt.nconf_in) then
4159           call make_ranvar(nconfr,m,idum)
4160           write (iout,*) nconfr-nconf_in,&
4161            " conformations have been generated randomly."
4162         endif
4163        else
4164         nconfr=nconfr*2
4165         call from_int(nconfr,m,idum)
4166 !       call from_pdb(nconfr,idum)
4167        endif
4168        write (iout,*) 'Exitted from make_ranvar nconfr=',nconfr
4169        write (*,*) 'Exitted from make_ranvar nconfr=',nconfr
4170        do m=1,nconfr
4171           write (iout,*) 'Initial conformation',m
4172           write(iout,'(8f10.4)') (rad2deg*dihang_in(1,j,1,m),j=2,nres-1)
4173           write(iout,'(8f10.4)') (rad2deg*dihang_in(2,j,1,m),j=2,nres-1)
4174           write(iout,'(8f10.4)') (rad2deg*dihang_in(3,j,1,m),j=2,nres-1)
4175           write(iout,'(8f10.4)') (rad2deg*dihang_in(4,j,1,m),j=2,nres-1)
4176        enddo 
4177        write(iout,*)'Calling FEEDIN NCONF',nconfr
4178        time1i=MPI_WTIME()
4179        call feedin(nconfr,nft)
4180        write (iout,*) ' Time for first bank min.',MPI_WTIME()-time1i
4181        call  history_append
4182         write(icsa_history,*) jlee,nft,nbank
4183         write(icsa_history,851) (etot(i),i=1,nconfr)
4184         write(icsa_history,850) (rmsn(i),i=1,nconfr)
4185         write(icsa_history,850) (pncn(i),i=1,nconfr)
4186         write(icsa_history,*)
4187        close(icsa_history)
4188       ELSE
4189 ! To minimize input conformation (bank conformation)
4190 ! Output to $mol.reminimized   
4191        if (irestart.lt.0) then 
4192         call reminimize(jlee)
4193         return
4194        endif
4195        if (irestart.eq.1) goto 111
4196 !  soldier - perform energy minimization
4197  334   call minim_jlee
4198       ENDIF
4199
4200 !cccccccccccccccccccccccccccccccccc
4201 ! need to syncronize all procs
4202       call mpi_barrier(CG_COMM,ierr)
4203       if (ierr.ne.0) then
4204        print *, ' cannot synchronize MPI'
4205        stop
4206       endif
4207 !cccccccccccccccccccccccccccccccccc
4208
4209       IF (ME.EQ.KING) THEN
4210
4211 !      print *,"ok after minim"
4212       nstep=nstep+nconf
4213       if(irestart.eq.2) then
4214        nbank=nbank+nconf
4215 !      ntbank=ntbank+nconf
4216        if(ntbank.gt.ntbankm) ntbank=ntbankm
4217       endif
4218 !      print *,"ok before indexx"
4219       if(iref.eq.0) then
4220        call indexx(nconfr,etot,indx)
4221       else
4222 ! cc/al 7/6/00
4223        do k=1,nconfr
4224          indx(k)=k
4225        enddo
4226        call indexx(nconfr-nconf_in,rmsn(nconf_in+1),indx(nconf_in+1))
4227        do k=nconf_in+1,nconfr
4228          indx(k)=indx(k)+nconf_in
4229        enddo
4230 ! cc/al
4231 !       call indexx(nconfr,rmsn,indx)
4232       endif
4233 !      print *,"ok after indexx"
4234       do im=1,nconf
4235        m=indx(im)
4236        if (m.gt.mxio .or. m.lt.1) then
4237          write (iout,*) 'Dimension ERROR in TOGEHER: IM',im,' M',m
4238          call mpi_abort(mpi_comm_world,ierror,ierrcode)
4239        endif
4240        jbank(im+nbank-nconf)=0
4241        bene(im+nbank-nconf)=etot(m)
4242        rene(im+nbank-nconf)=etot(m)
4243 !!bankt       btene(im)=etot(m)
4244 !
4245        brmsn(im+nbank-nconf)=rmsn(m)
4246        bpncn(im+nbank-nconf)=pncn(m)
4247        rrmsn(im+nbank-nconf)=rmsn(m)
4248        rpncn(im+nbank-nconf)=pncn(m)
4249        if (im+nbank-nconf.gt.mxio .or. im+nbank-nconf.lt.1) then
4250          write (iout,*) 'Dimension ERROR in TOGEHER: IM',im,&
4251          ' NBANK',nbank,' NCONF',nconf,' IM+NBANK-NCONF',im+nbank-nconf
4252          call mpi_abort(mpi_comm_world,ierror,ierrcode)
4253        endif
4254        do k=1,numch
4255         do j=2,nres-1
4256          do i=1,4
4257           bvar(i,j,k,im+nbank-nconf)=dihang(i,j,k,m)
4258           rvar(i,j,k,im+nbank-nconf)=dihang(i,j,k,m)
4259 !!bankt          btvar(i,j,k,im)=dihang(i,j,k,m)
4260 !
4261          enddo
4262         enddo
4263        enddo
4264        if(iref.eq.1) then
4265         if(brmsn(im+nbank-nconf).gt.rmscut.or. &
4266            bpncn(im+nbank-nconf).lt.pnccut) ibank(im+nbank-nconf)=9
4267        endif
4268        if(vdisulf) then
4269            bvar_ns(im+nbank-nconf)=ns-2*nss
4270            k=0
4271            do i=1,ns
4272              j=1
4273              do while( iss(i).ne.ihpb(j)-nres .and. &
4274                        iss(i).ne.jhpb(j)-nres .and. j.le.nss)
4275               j=j+1 
4276              enddo
4277              if (j.gt.nss) then            
4278                k=k+1
4279                bvar_s(k,im+nbank-nconf)=iss(i)
4280              endif
4281            enddo
4282        endif
4283        bvar_nss(im+nbank-nconf)=nss
4284        do i=1,nss
4285            bvar_ss(1,i,im+nbank-nconf)=ihpb(i)
4286            bvar_ss(2,i,im+nbank-nconf)=jhpb(i)
4287        enddo
4288       enddo
4289       ENDIF
4290
4291   111 continue
4292
4293       IF (ME.EQ.KING) THEN
4294
4295       call find_max
4296       call find_min
4297
4298       call get_diff
4299       if(nbank.eq.nconf.and.irestart.eq.0) then
4300        adif=avedif
4301       endif
4302
4303       cutdif=adif/cut1
4304       ctdif1=adif/cut2
4305
4306 !d      print *,"adif,xctdif,cutdifr"
4307 !d      print *,adif,xctdif,cutdifr
4308        nst=ntotal/ntrial/nseed
4309        xctdif=(cutdif/ctdif1)**(-1.0/nst)
4310        if(irestart.ge.1) call estimate_cutdif(adif,xctdif,cutdifr)
4311 !       print *,"ok after estimate"
4312
4313       irestart=0
4314
4315        call write_rbank(jlee,adif,nft)
4316        call write_bank(jlee,nft)
4317 !!bankt       call write_bankt(jlee,nft)
4318 !       call write_bank1(jlee)
4319        call  history_append
4320         write(icsa_history,*) "xctdif: ", xctdif,nst,adif/cut1,ctdif1
4321         write(icsa_history,851) (bene(i),i=1,nbank)
4322         write(icsa_history,850) (brmsn(i),i=1,nbank)
4323         write(icsa_history,850) (bpncn(i),i=1,nbank)
4324         close(icsa_history)
4325   850 format(10f8.3)
4326   851 format(5e15.6)
4327
4328       ifar=nseed/4*3+1
4329       ifar=nseed+1
4330       ENDIF
4331     
4332
4333       finished=.false.
4334       iter = 0
4335       irecv = 0
4336       isent =0
4337       ifrom= 0
4338       time0i=MPI_WTIME()
4339       time1i=time0i
4340       time_start_c=time0i
4341       if (.not.sync_iter) then 
4342         time_start=time0i
4343         nft00=nft
4344       else
4345         sync_iter=.false.
4346       endif
4347       nft00_c=nft
4348       nft0i=nft
4349 !cccccccccccccccccccccccccccccccccccccc
4350       do while (.not. finished)
4351 !cccccccccccccccccccccccccccccccccccccc
4352 !rc      print *,"iter ", iter,' isent=',isent
4353
4354       IF (ME.EQ.KING) THEN
4355 !  start energy minimization
4356
4357        if (isent.eq.0) then
4358 !  king-emperor - select seeds & make var & feed input
4359 !d        print *,'generating new conf',ntrial,MPI_WTIME()
4360         call select_is(nseed,ifar,idum)
4361
4362         open(icsa_seed,file=csa_seed,status="old")
4363         write(icsa_seed,39)  &
4364           jlee,icycle,nstep,(is(i),bene(is(i)),i=1,nseed)
4365         close(icsa_seed)
4366         call  history_append
4367         write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,&
4368          ebmin,ebmax,nft,iuse,nbank,ntbank
4369         close(icsa_history)
4370
4371          
4372
4373         call make_var(ntry,idum,iter)
4374 !d        print *,'new trial generated',ntrial,MPI_WTIME()
4375            time2i=MPI_WTIME()
4376            write (iout,'(a20,i4,f12.2)') &
4377              'Time for make trial',iter+1,time2i-time1i
4378        endif
4379
4380 !rc        write(iout,*)'1:Calling FEEDIN NTRY',NTRY,' ntrial',ntrial
4381 !rc        call feedin(ntry,nft)
4382
4383        isent=isent+1
4384        if (isent.ge.nodes.or.iter.gt.0)  then
4385 !t            print *,'waiting ',MPI_WTIME()
4386             irecv=irecv+1
4387             call recv(0,ifrom,xout,eout,ind,timeout)
4388 !t            print *,'   ',irecv,' received from',ifrom,MPI_WTIME()
4389        else
4390             ifrom=ifrom+1
4391        endif
4392
4393 !t            print *,'sending to',ifrom,MPI_WTIME()
4394        call send(isent,ifrom,iter)
4395 !t            print *,isent,' sent ',MPI_WTIME()
4396
4397 ! store results -----------------------------------------------
4398        if (isent.ge.nodes.or.iter.gt.0)  then
4399          nft=nft+ind(3)
4400          movernx(irecv)=iabs(ind(5))
4401          call getx(ind,xout,eout,cout,rad,iw_pdb,irecv)
4402          if(vdisulf) then
4403              nss_out(irecv)=nss
4404              do i=1,nss
4405                iss_out(i,irecv)=ihpb(i)
4406                jss_out(i,irecv)=jhpb(i)  
4407              enddo
4408          endif
4409          if(iw_pdb.gt.0) &
4410                 call write_csa_pdb(xout,eout,nft,irecv,iw_pdb)
4411        endif
4412 !--------------------------------------------------------------
4413        if (isent.eq.ntry) then
4414            time1i=MPI_WTIME()
4415            write (iout,'(a18,f12.2,a14,f10.2)') &
4416              'Nonsetup time     ',time1i-time_start_c,&
4417              ' sec, Eval/s =',(nft-nft00_c)/(time1i-time_start_c)
4418            write (iout,'(a14,i4,f12.2,a14,f10.2)') &
4419              'Time for iter ',iter+1,time1i-time0i,&
4420              ' sec, Eval/s =',(nft-nft0i)/(time1i-time0i)
4421            time0i=time1i
4422            nft0i=nft
4423            cutdif=cutdif*xctdif
4424            if(cutdif.lt.ctdif1) cutdif=ctdif1
4425            if (iter.eq.0) then
4426               print *,'UPDATING ',ntry-nodes+1,irecv
4427               write(iout,*) 'UPDATING ',ntry-nodes+1
4428               iter=iter+1
4429 !----------------- call update(ntry-nodes+1) -------------------
4430               nstep=nstep+ntry-nseed-(nodes-1)
4431               call refresh_bank(ntry-nodes+1)
4432 !!bankt              call refresh_bankt(ntry-nodes+1)
4433            else
4434 !----------------- call update(ntry) ---------------------------
4435               iter=iter+1
4436               print *,'UPDATING ',ntry,irecv
4437               write(iout,*) 'UPDATING ',ntry
4438               nstep=nstep+ntry-nseed
4439               call refresh_bank(ntry)
4440 !!bankt              call refresh_bankt(ntry)
4441            endif         
4442 !----------------------------------------------------------------- 
4443
4444            call write_bank(jlee,nft)
4445 !!bankt           call write_bankt(jlee,nft)
4446            call find_min
4447
4448            time1i=MPI_WTIME()
4449            write (iout,'(a20,i4,f12.2)') &
4450              'Time for refresh ',iter,time1i-time0i
4451
4452            if(ebmin.lt.estop) finished=.true.
4453            if(icycle.gt.icmax) then
4454                call write_bank1(jlee)
4455                do i=1,nbank
4456                  ibank(i)=2
4457                  ibank(i)=1
4458                enddo
4459                nbank=nbank+nconf
4460                if(nbank.gt.1000) then 
4461                    finished=.true.
4462                else
4463 !rc                   goto 333
4464                    sync_iter=.true.
4465                endif
4466            endif
4467            if(nstep.gt.nstmax) finished=.true.
4468
4469            if(finished.or.sync_iter) then
4470             do ij=1,nodes-1
4471               call recv(1,ifrom,xout,eout,ind,timeout)
4472               if (timeout) then
4473                 nstep=nstep+ij-1
4474                 print *,'ERROR worker is not responding'
4475                 write(iout,*) 'ERROR worker is not responding'
4476                 time1i=MPI_WTIME()-time_start_c
4477                 print *,'End of cycle, master time for ',iter,' iters ',&
4478                    time1i,'sec, Eval/s ',(nft-nft00_c)/time1i
4479                 write (iout,*) 'End of cycle, master time for ',iter,&
4480                    ' iters ',time1i,' sec'
4481                 write (iout,*) 'Total eval/s ',(nft-nft00_c)/time1i
4482                 print *,'UPDATING ',ij-1
4483                 write(iout,*) 'UPDATING ',ij-1
4484                 call flush(iout)
4485                 call refresh_bank(ij-1)
4486 !!bankt                call refresh_bankt(ij-1)
4487                 goto 1002
4488               endif
4489 !              print *,'node ',ifrom,' finished ',ij,nft
4490               write(iout,*) 'node ',ifrom,' finished ',ij,nft
4491               call flush(iout)
4492               nft=nft+ind(3)
4493               movernx(ij)=iabs(ind(5))
4494               call getx(ind,xout,eout,cout,rad,iw_pdb,ij)
4495               if(vdisulf) then
4496                nss_out(ij)=nss
4497                do i=1,nss
4498                  iss_out(i,ij)=ihpb(i)
4499                  jss_out(i,ij)=jhpb(i)  
4500                enddo
4501               endif
4502               if(iw_pdb.gt.0) &
4503                 call write_csa_pdb(xout,eout,nft,ij,iw_pdb)
4504             enddo
4505             nstep=nstep+nodes-1
4506 !rc            print *,'---------bcast finished--------',finished
4507             time1i=MPI_WTIME()-time_start_c
4508             print *,'End of cycle, master time for ',iter,' iters ',&
4509                    time1i,'sec, Eval/s ',(nft-nft00_c)/time1i
4510             write (iout,*) 'End of cycle, master time for ',iter,&
4511                    ' iters ',time1i,' sec'
4512             write (iout,*) 'Total eval/s ',(nft-nft00_c)/time1i
4513
4514 !timeout            call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr)
4515 !timeout            call mpi_bcast(sync_iter,1,mpi_logical,0,
4516 !timeout     &                                              CG_COMM,ierr)
4517             do ij=1,nodes-1 
4518                tstart=MPI_WTIME()
4519                call mpi_issend(finished,1,mpi_logical,ij,idchar,&
4520                    CG_COMM,ireq,ierr)
4521                call mpi_issend(sync_iter,1,mpi_logical,ij,idchar,&
4522                    CG_COMM,ireq2,ierr)
4523                flag=.false.  
4524                timeout1=.false.
4525                do while(.not. (flag .or. timeout1))
4526                  call MPI_TEST(ireq2,flag,muster,ierr)
4527                  tend1=MPI_WTIME()
4528                  if(tend1-tstart.gt.60) then 
4529                   print *,'ERROR worker ',ij,' is not responding'
4530                   write(iout,*) 'ERROR worker ',ij,' is not responding'
4531                   timeout1=.true.
4532                  endif
4533                enddo
4534                if(timeout1) then
4535                 write(iout,*) 'worker ',ij,' NOT OK ',tend1-tstart
4536                 timeout=.true.
4537                else
4538                 write(iout,*) 'worker ',ij,' OK ',tend1-tstart
4539                endif
4540             enddo
4541             print *,'UPDATING ',nodes-1,ij
4542             write(iout,*) 'UPDATING ',nodes-1
4543             call refresh_bank(nodes-1)
4544 !!bankt            call refresh_bankt(nodes-1)
4545  1002       continue
4546             call write_bank(jlee,nft)
4547 !!bankt            call write_bankt(jlee,nft)
4548             call find_min
4549
4550             do i=0,mxmv  
4551               do j=1,3
4552                 nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j)
4553                 nstatnx(i,j)=0
4554               enddo
4555             enddo
4556
4557             write(iout,*)'### Total stats:'
4558             do i=0,mxmv
4559              if(nstatnx_tot(i,1).ne.0) then
4560               if (i.le.9) then
4561               write(iout,'(a5,i1,a7,i6,a7,i4,a5,i4,a5,f5.1)') &
4562               '### N',i,' total=',nstatnx_tot(i,1),&
4563             ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),'%acc',&
4564              (nstatnx_tot(i,2)+nstatnx_tot(i,3))*100.0/nstatnx_tot(i,1)
4565               else
4566               write(iout,'(a4,i2,a7,i6,a7,i4,a5,i4,a5,f5.1)') &
4567               '###N',i,' total=',nstatnx_tot(i,1),&
4568             ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),'%acc',&
4569              (nstatnx_tot(i,2)+nstatnx_tot(i,3))*100.0/nstatnx_tot(i,1)
4570               endif
4571              else
4572               if (i.le.9) then
4573               write(iout,'(a5,i1,a7,i6,a7,i4,a5,i4,a5,f5.1)') &
4574                 '### N',i,' total=',nstatnx_tot(i,1),&
4575                 ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),&
4576                 ' %acc',0.0
4577               else
4578               write(iout,'(a4,i2,a7,i6,a7,i4,a5,i4,a5,f5.1)') &
4579                 '###N',i,' total=',nstatnx_tot(i,1),&
4580                 ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),&
4581                 ' %acc',0.0
4582               endif
4583              endif
4584             enddo
4585
4586            endif
4587            if(sync_iter) goto 331
4588
4589    39 format(2i3,i7,5(i4,f8.3,1x),19(/,13x,5(i4,f8.3,1x)))
4590    40 format(2i2,i8,f8.1,2i4,2(1pe14.5),i10,3i4)
4591    43 format(10i8)
4592    44 format('jlee =',i3,':',4f10.1,' E =',f8.3,i7,i10)
4593
4594            isent=0
4595            irecv=0
4596        endif
4597       ELSE
4598 !  soldier - perform energy minimization
4599         call minim_jlee
4600         print *,'End of minim, proc',me,'time ',MPI_WTIME()-time_start
4601         write (iout,*) 'End of minim, proc',me,'time ',&
4602                         MPI_WTIME()-time_start
4603         call flush(iout)
4604 !timeout        call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr)
4605 !timeout        call mpi_bcast(sync_iter,1,mpi_logical,0,CG_COMM,ierr)
4606          call mpi_recv(finished,1,mpi_logical,0,idchar,&
4607                        CG_COMM,muster,ierr)             
4608          call mpi_recv(sync_iter,1,mpi_logical,0,idchar,&
4609                        CG_COMM,muster,ierr)             
4610         if(sync_iter) goto 331
4611       ENDIF
4612
4613 !cccccccccccccccccccccccccccccccccccccc
4614       enddo
4615 !cccccccccccccccccccccccccccccccccccccc
4616
4617       IF (ME.EQ.KING) THEN
4618         call  history_append
4619         write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,&
4620         ebmin,ebmax,nft,iuse,nbank,ntbank
4621
4622         write(icsa_history,44) jlee,0.0,0.0,0.0,&
4623          0.0,ebmin,nstep,nft
4624         write(icsa_history,*)
4625        close(icsa_history)
4626
4627        time1i=MPI_WTIME()-time_start
4628        print *,'End of RUN, master time ',&
4629                    time1i,'sec, Eval/s ',(nft-nft00)/time1i
4630        write (iout,*) 'End of RUN, master time  ',&
4631                    time1i,' sec'
4632        write (iout,*) 'Total eval/s ',(nft-nft00)/time1i
4633
4634        if(timeout) then 
4635         write(iout,*) '!!!! ERROR worker was not responding'
4636         write(iout,*) '!!!! cannot finish work normally'
4637         write(iout,*) 'Processor0 is calling MPI_ABORT'
4638         print *,'!!!! ERROR worker was not responding'
4639         print *,'!!!! cannot finish work normally'
4640         print *,'Processor0 is calling MPI_ABORT'
4641         call flush(iout)
4642         call mpi_abort(mpi_comm_world, 111, ierr)
4643        endif
4644       ENDIF
4645
4646 !ccccccccccccccccccccccccccccc
4647   300 continue
4648 !ccccccccccccccccccccccccccccc
4649
4650       return
4651       end subroutine together
4652 !-----------------------------------------------------------------------------
4653       subroutine feedin(nconf,nft)
4654
4655       use MPI_data
4656       use geometry_data, only:nvar
4657       use io, only:write_csa_pdb
4658 !  sends out starting conformations and receives results of energy minimization
4659 !      implicit real*8 (a-h,o-z)
4660 !      include 'DIMENSIONS'
4661 !      include 'COMMON.VAR'
4662 !      include 'COMMON.IOUNITS'
4663 !      include 'COMMON.CONTROL'
4664       include 'mpif.h'
4665       real(kind=8),dimension(6*nres) :: xin,xout        !(maxvar) (maxvar=6*maxres)
4666       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
4667       real(kind=8),dimension(2) :: cout
4668       integer,dimension(9) :: ind
4669       integer,dimension(12) :: info
4670       integer,dimension(mpi_status_size) :: muster
4671 !      include 'COMMON.SETUP'
4672       real(kind=8),parameter :: rad=1.745329252d-2
4673       integer :: j,nconf,nft,mm,n,ierror,ierrcode,ierr,iw_pdb,&
4674             man
4675
4676       print *,'FEEDIN: NCONF=',nconf
4677       mm=0
4678 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4679       if (nconf .lt. nodes-1) then
4680         write (*,*) 'FATAL ERROR in FEEDIN, nconf < nodes -1',&
4681          nconf,nodes-1 
4682         write (iout,*) 'FATAL ERROR in FEEDIN, nconf < nodes -1',&
4683          nconf,nodes-1 
4684         call mpi_abort(mpi_comm_world,ierror,ierrcode)
4685       endif
4686       do n=1,nconf
4687 !  pull out external and internal variables for next start
4688         call putx(xin,n,rad)
4689 !        write (iout,*) 'XIN from FEEDIN N=',n
4690 !        write(iout,'(8f10.4)') (xin(j),j=1,nvar)
4691         mm=mm+1
4692         if (mm.lt.nodes) then
4693 !  feed task to soldier
4694 !       print *, ' sending input for start # ',n
4695          info(1)=n
4696          info(2)=-1
4697          info(3)=0
4698          info(4)=0
4699          info(5)=0
4700          info(6)=0
4701          call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,&
4702                         ierr)
4703          call mpi_send(xin,nvar,mpi_double_precision,mm,&
4704                         idreal,CG_COMM,ierr)
4705         else
4706 !  find an available soldier
4707          call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,&
4708                        CG_COMM,muster,ierr)
4709 !        print *, ' receiving output from start # ',ind(1)
4710          man=muster(mpi_source)
4711 !  receive final energies and variables
4712          nft=nft+ind(3)
4713          call mpi_recv(eout,1,mpi_double_precision,&
4714                        man,idreal,CG_COMM,muster,ierr)
4715 !         print *,eout 
4716 #ifdef CO_BIAS
4717          call mpi_recv(co,1,mpi_double_precision,&
4718                        man,idreal,CG_COMM,muster,ierr)
4719          write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
4720 #endif
4721          call mpi_recv(xout,nvar,mpi_double_precision,&
4722                        man,idreal,CG_COMM,muster,ierr)
4723 !         print *,nvar , ierr
4724 !  feed next task to soldier
4725 !        print *, ' sending input for start # ',n
4726          info(1)=n
4727          info(2)=-1
4728          info(3)=0
4729          info(4)=0
4730          info(5)=0
4731          info(6)=0
4732          info(7)=0
4733          info(8)=0
4734          info(9)=0
4735          call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,&
4736                         ierr)
4737          call mpi_send(xin,nvar,mpi_double_precision,man,&
4738                         idreal,CG_COMM,ierr)
4739 !  retrieve latest results
4740          call getx(ind,xout,eout,cout,rad,iw_pdb,ind(1))
4741          if(iw_pdb.gt.0) &
4742               call write_csa_pdb(xout,eout,nft,ind(1),iw_pdb)
4743         endif
4744       enddo
4745 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4746 !  no more input
4747 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4748       do j=1,nodes-1
4749 !  wait for a soldier
4750        call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,&
4751                      CG_COMM,muster,ierr)
4752 !rc       if (ierr.ne.0) go to 30
4753 !      print *, ' receiving output from start # ',ind(1)
4754        man=muster(mpi_source)
4755 !  receive final energies and variables
4756        nft=nft+ind(3)
4757        call mpi_recv(eout,1,&
4758                      mpi_double_precision,man,idreal,&
4759                      CG_COMM,muster,ierr)
4760 !       print *,eout
4761 #ifdef CO_BIAS
4762          call mpi_recv(co,1,mpi_double_precision,&
4763                        man,idreal,CG_COMM,muster,ierr)
4764          write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
4765 #endif
4766 !rc       if (ierr.ne.0) go to 30
4767        call mpi_recv(xout,nvar,mpi_double_precision,&
4768                      man,idreal,CG_COMM,muster,ierr)
4769 !       print *,nvar , ierr
4770 !rc       if (ierr.ne.0) go to 30
4771 !  halt soldier
4772        info(1)=0
4773        info(2)=-1
4774        info(3)=0 
4775        info(4)=0
4776        info(5)=0
4777        info(6)=0
4778        call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,&
4779                       ierr)
4780 !  retrieve results
4781        call getx(ind,xout,eout,cout,rad,iw_pdb,ind(1))
4782        if(iw_pdb.gt.0) &
4783                 call write_csa_pdb(xout,eout,nft,ind(1),iw_pdb)
4784       enddo
4785 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4786       return
4787    10 print *, ' dispatching error'
4788       call mpi_abort(mpi_comm_world,ierror,ierrcode)
4789       return
4790    20 print *, ' communication error'
4791       call mpi_abort(mpi_comm_world,ierror,ierrcode)
4792       return
4793    30 print *, ' receiving error'
4794       call mpi_abort(mpi_comm_world,ierror,ierrcode)
4795
4796       return
4797       end subroutine feedin
4798 !-----------------------------------------------------------------------------
4799       subroutine getx(ind,xout,eout,cout,rad,iw_pdb,k)
4800
4801       use geometry_data
4802       use energy_data
4803       use compare, only: contact_fract
4804       use MPI_data
4805       include 'mpif.h'
4806 !  receives and stores data from soldiers
4807 !      implicit real*8 (a-h,o-z)
4808 !      include 'DIMENSIONS'
4809 !      include 'COMMON.IOUNITS'
4810 !      include 'COMMON.CSA'
4811 !      include 'COMMON.BANK'
4812 !      include 'COMMON.VAR'
4813 !      include 'COMMON.CHAIN'
4814 !      include 'COMMON.CONTACTS'
4815       integer,dimension(9) :: ind
4816       real(kind=8),dimension(6*nres) :: xout    !(maxvar) (maxvar=6*maxres)
4817       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
4818 !jlee
4819       real(kind=8) :: przes(3),obr(3,3),cout(2)
4820       logical :: non_conv
4821       integer :: iw_pdb,k,j,ierror,ierrcode
4822       real(kind=8) :: rad,co
4823 !jlee
4824       iw_pdb=2
4825       if (k.gt.mxio .or. k.lt.1) then 
4826         write (iout,*) &
4827          'ERROR - dimensions of ANGMIN have been exceeded K=',k
4828         call mpi_abort(mpi_comm_world,ierror,ierrcode)
4829       endif
4830 !  store ind()
4831       do j=1,9
4832        indb(k,j)=ind(j)
4833       enddo
4834 !  store energies
4835       etot(k)=eout(1)
4836 !  retrieve dihedral angles etc
4837       call var_to_geom(nvar,xout)
4838       do j=2,nres-1
4839         dihang(1,j,1,k)=theta(j+1)
4840         dihang(2,j,1,k)=phi(j+2)
4841         dihang(3,j,1,k)=alph(j)
4842         dihang(4,j,1,k)=omeg(j)
4843       enddo
4844       dihang(2,nres-1,1,k)=0.0d0
4845 !jlee
4846       if(iref.eq.0) then 
4847        iw_pdb=1
4848 !d       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4)') 
4849 !d     &      ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' mv ',
4850 !d     &      ind(5),ind(4)
4851        return
4852       endif
4853       call chainbuild
4854 !     call dihang_to_c(dihang(1,1,1,k))
4855 !     call fitsq(rms,c(1,1),crefjlee(1,1),nres,przes,obr,non_conv)
4856 !     call fitsq(rms,c(1,2),crefjlee(1,2),nres-1,przes,obr,non_conv)
4857 !           call fitsq(rms,c(1,nstart_seq),crefjlee(1,nstart_sup),
4858 !    &                 nsup,przes,obr,non_conv)
4859 !      rmsn(k)=dsqrt(rms)
4860
4861        call rmsd_csa(rmsn(k))
4862        call contact(.false.,ncont,icont,co)
4863        pncn(k)=contact_fract(ncont,ncont_ref,icont,icont_ref)     
4864
4865 !d       write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a5
4866 !d     &      ,0pf5.2,a5,f5.1,a,f6.3,a4,i3,i4)') 
4867 !d     &    ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' rms ',
4868 !d     &    rmsn(k),' %NC ',pncn(k)*100,' cont.order',co,' mv ',
4869 !d     &    ind(5),ind(4)
4870
4871  
4872       if (rmsn(k).gt.rmscut.or.pncn(k).lt.pnccut) iw_pdb=0
4873       return
4874       end subroutine getx
4875 !-----------------------------------------------------------------------------
4876       subroutine putx(xin,n,rad)
4877
4878       use geometry_data
4879 !  gets starting variables
4880 !      implicit real*8 (a-h,o-z)
4881 !      include 'DIMENSIONS'
4882 !      include 'COMMON.CSA'
4883 !      include 'COMMON.BANK'
4884 !      include 'COMMON.VAR'
4885 !      include 'COMMON.CHAIN'
4886 !      include 'COMMON.IOUNITS'
4887       integer :: n,m,j
4888       real(kind=8),dimension(6*nres) :: xin     !(maxvar) (maxvar=6*maxres)
4889       real(kind=8) :: rad
4890
4891 !  pull out starting values for variables
4892 !       write (iout,*)'PUTX: N=',n
4893       do m=1,numch
4894 !        write (iout,'(8f10.4)') (dihang_in(1,j,m,n),j=2,nres-1)
4895 !        write (iout,'(8f10.4)') (dihang_in(2,j,m,n),j=2,nres-1)
4896 !        write (iout,'(8f10.4)') (dihang_in(3,j,m,n),j=2,nres-1)
4897 !        write (iout,'(8f10.4)') (dihang_in(4,j,m,n),j=2,nres-1)
4898        do j=2,nres-1
4899         theta(j+1)=dihang_in(1,j,m,n)
4900         phi(j+2)=dihang_in(2,j,m,n)
4901         alph(j)=dihang_in(3,j,m,n)
4902         omeg(j)=dihang_in(4,j,m,n)
4903        enddo
4904       enddo
4905 !  set up array of variables
4906       call geom_to_var(nvar,xin)
4907 !       write (iout,*) 'xin in PUTX N=',n 
4908 !       call intout
4909 !       write (iout,'(8f10.4)') (xin(i),i=1,nvar) 
4910       return
4911       end subroutine putx
4912 !-----------------------------------------------------------------------------
4913       subroutine putx2(xin,iff,n)
4914
4915       use geometry_data
4916 !  gets starting variables
4917 !      implicit real*8 (a-h,o-z)
4918 !      include 'DIMENSIONS'
4919 !      include 'COMMON.CSA'
4920 !      include 'COMMON.BANK'
4921 !      include 'COMMON.VAR'
4922 !      include 'COMMON.CHAIN'
4923 !      include 'COMMON.IOUNITS'
4924       integer :: n,m,j,i
4925       real(kind=8),dimension(6*nres) :: xin     !(maxvar) (maxvar=6*maxres)
4926       integer,dimension(nres) :: iff    !(maxres)
4927
4928 !  pull out starting values for variables
4929       do m=1,numch
4930        do j=2,nres-1
4931         theta(j+1)=dihang_in2(1,j,m,n)
4932         phi(j+2)=dihang_in2(2,j,m,n)
4933         alph(j)=dihang_in2(3,j,m,n)
4934         omeg(j)=dihang_in2(4,j,m,n)
4935        enddo
4936       enddo
4937 !  set up array of variables
4938       call geom_to_var(nvar,xin)
4939        
4940       do i=1,nres
4941         iff(i)=iff_in(i,n)
4942       enddo
4943       return
4944       end subroutine putx2
4945 !-----------------------------------------------------------------------------
4946       subroutine prune_bank(p_cut)
4947
4948       use MPI_data
4949 !      implicit real*8 (a-h,o-z)
4950 !      include 'DIMENSIONS'
4951       include 'mpif.h'
4952 !      include 'COMMON.CSA'
4953 !      include 'COMMON.BANK'
4954 !      include 'COMMON.IOUNITS'
4955 !      include 'COMMON.CHAIN'
4956 !      include 'COMMON.TIME1'
4957 !      include 'COMMON.SETUP'
4958       integer :: k,j,i,m,ip,nprune
4959       real(kind=8) :: p_cut,diff,ddmin
4960 !---------------------------
4961 ! This subroutine prunes bank conformations using p_cut
4962 !---------------------------
4963
4964       nprune=0
4965       nprune=nprune+1
4966       m=1
4967       do k=1,numch
4968        do j=2,nres-1
4969         do i=1,4
4970          dihang(i,j,k,nprune)=bvar(i,j,k,m)
4971         enddo
4972        enddo
4973       enddo
4974       bene(nprune)=bene(m)
4975       brmsn(nprune)=brmsn(m)
4976       bpncn(nprune)=bpncn(m) 
4977
4978       do m=2,nbank
4979        ddmin=9.d190
4980        do ip=1,nprune
4981         call get_diff12(dihang(1,1,1,ip),bvar(1,1,1,m),diff) 
4982         if(diff.lt.p_cut) goto 100
4983         if(diff.lt.ddmin) ddmin=diff
4984        enddo
4985        nprune=nprune+1
4986        do k=1,numch
4987         do j=2,nres-1
4988          do i=1,4
4989           dihang(i,j,k,nprune)=bvar(i,j,k,m)
4990          enddo
4991         enddo
4992        enddo
4993        bene(nprune)=bene(m)
4994        brmsn(nprune)=brmsn(m)
4995        bpncn(nprune)=bpncn(m)
4996   100  continue
4997        write (iout,*) 'Pruning :',m,nprune,p_cut,ddmin
4998       enddo
4999       nbank=nprune
5000       print *, 'Pruning :',m,nprune,p_cut
5001       call write_bank(0,0)
5002
5003       return
5004       end subroutine prune_bank
5005 !-----------------------------------------------------------------------------
5006       subroutine reminimize(jlee)
5007
5008       use MPI_data
5009 !      implicit real*8 (a-h,o-z)
5010 !      include 'DIMENSIONS'
5011       include 'mpif.h'
5012 !      include 'COMMON.CSA'
5013 !      include 'COMMON.BANK'
5014 !      include 'COMMON.IOUNITS'
5015 !      include 'COMMON.CHAIN'
5016 !      include 'COMMON.TIME1'
5017 !      include 'COMMON.SETUP'
5018       integer :: i,j,k,jlee,index,nft,ntry
5019 !---------------------------
5020 ! This subroutine re-minimizes bank conformations:
5021 !---------------------------
5022
5023        ntry=nbank
5024
5025        call find_max
5026        call find_min
5027
5028        if (me.eq.king) then
5029         open(icsa_history,file=csa_history,status="old")
5030          write(icsa_history,*) "Re-minimization",nodes,"nodes"
5031          write(icsa_history,851) (bene(i),i=1,nbank)
5032          write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,&
5033           ebmin,ebmax,nft,iuse,nbank,ntbank
5034         close(icsa_history)
5035         do index=1,ntry
5036          do k=1,numch
5037           do j=2,nres-1
5038            do i=1,4
5039             dihang_in(i,j,k,index)=bvar(i,j,k,index)
5040            enddo
5041           enddo
5042          enddo
5043         enddo
5044         nft=0
5045         call feedin(ntry,nft)
5046        else
5047         call minim_jlee
5048        endif
5049
5050        call find_max
5051        call find_min
5052
5053        if (me.eq.king) then
5054         do i=1,ntry
5055          call replace_bvar(i,i)
5056         enddo
5057         open(icsa_history,file=csa_history,status="old")
5058          write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,&
5059           ebmin,ebmax,nft,iuse,nbank,ntbank
5060          write(icsa_history,851) (bene(i),i=1,nbank)
5061         close(icsa_history)
5062         call write_bank_reminimized(jlee,nft)
5063        endif
5064
5065    40 format(2i2,i8,f8.1,2i4,2(1pe14.5),i10,3i4)
5066   851 format(5e15.6)
5067   850 format(5e15.10)
5068 !  850 format(10f8.3)
5069
5070       return
5071       end subroutine reminimize
5072 !-----------------------------------------------------------------------------
5073       subroutine send(n,mm,it)
5074
5075       use MPI_data
5076       use geometry_data, only: nvar
5077       use control_data, only: vdisulf
5078 !  sends out starting conformation for minimization
5079 !      implicit real*8 (a-h,o-z)
5080 !      include 'DIMENSIONS'
5081 !      include 'COMMON.VAR'
5082 !      include 'COMMON.IOUNITS'
5083 !      include 'COMMON.CONTROL'
5084 !      include 'COMMON.BANK'
5085 !      include 'COMMON.CHAIN'
5086       include 'mpif.h'
5087       real(kind=8),dimension(6*nres) :: xin,xout,xin2   !(maxvar) (maxvar=6*maxres)
5088       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
5089       real(kind=8),dimension(2) :: cout
5090       integer,dimension(9) :: ind
5091       integer,dimension(nres) :: iff    !(maxres)
5092       integer,dimension(12) :: info
5093       integer,dimension(mpi_status_size) :: muster
5094 !      include 'COMMON.SETUP'
5095       real(kind=8),parameter :: rad=1.745329252d-2
5096       integer :: n,mm,it,ierr
5097
5098       if (isend2(n).eq.0) then
5099 !  pull out external and internal variables for next start
5100         call putx(xin,n,rad)
5101         info(1)=n
5102         info(2)=it
5103         info(3)=movenx(n)
5104         info(4)=nss_in(n)
5105         info(5)=parent(1,n)
5106         info(6)=parent(2,n)
5107
5108         if (movenx(n).eq.14.or.movenx(n).eq.17) then
5109           info(7)=idata(1,n)
5110           info(8)=idata(2,n)
5111         else if (movenx(n).eq.16) then
5112           info(7)=idata(1,n)
5113           info(8)=idata(2,n)
5114           info(10)=idata(3,n)
5115           info(11)=idata(4,n)
5116           info(12)=idata(5,n)
5117         else
5118          info(7)=0
5119          info(8)=0
5120          info(10)=0
5121          info(11)=0
5122          info(12)=0
5123         endif
5124
5125         if (movenx(n).eq.15) then
5126          info(9)=parent(3,n)
5127         else
5128          info(9)=0
5129         endif
5130         call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,&
5131                         ierr)
5132         call mpi_send(xin,nvar,mpi_double_precision,mm,&
5133                         idreal,CG_COMM,ierr)
5134       else
5135 !  distfit & minimization for n7 move
5136         info(1)=-n
5137         info(2)=it
5138         info(3)=movenx(n)
5139         info(4)=nss_in(n)
5140         info(5)=parent(1,n)
5141         info(6)=parent(2,n)
5142         info(7)=0
5143         info(8)=0
5144         info(9)=0
5145         call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,&
5146                         ierr)
5147         call putx2(xin,iff,isend2(n))
5148         call mpi_send(xin,nvar,mpi_double_precision,mm,&
5149                         idreal,CG_COMM,ierr)
5150         call mpi_send(iff,nres,mpi_integer,mm,&
5151                         idint,CG_COMM,ierr)
5152         call putx(xin2,n,rad)
5153         call mpi_send(xin2,nvar,mpi_double_precision,mm,&
5154                         idreal,CG_COMM,ierr)
5155       endif 
5156       if (vdisulf.and.nss_in(n).ne.0) then
5157         call mpi_send(iss_in(1,n),nss_in(n),mpi_integer,mm,&
5158                         idint,CG_COMM,ierr)
5159         call mpi_send(jss_in(1,n),nss_in(n),mpi_integer,mm,&
5160                         idint,CG_COMM,ierr)
5161       endif
5162       return
5163       end subroutine send
5164 !-----------------------------------------------------------------------------
5165       subroutine recv(ihalt,man,xout,eout,ind,tout)
5166
5167       use MPI_data
5168       use energy_data
5169       use geometry_data, only: nvar
5170       use control_data, only: vdisulf
5171 !  receives results of energy minimization
5172 !      implicit real*8 (a-h,o-z)
5173 !      include 'DIMENSIONS'
5174 !      include 'COMMON.VAR'
5175 !      include 'COMMON.IOUNITS'
5176 !      include 'COMMON.CONTROL'
5177 !      include 'COMMON.SBRIDGE'
5178 !      include 'COMMON.BANK'
5179 !      include 'COMMON.CHAIN'
5180       include 'mpif.h'
5181       real(kind=8),dimension(6*nres) :: xin,xout        !(maxvar) (maxvar=6*maxres)
5182       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout
5183       real(kind=8),dimension(2) :: cout
5184       integer,dimension(9) :: ind
5185       integer,dimension(12) :: info
5186       integer,dimension(mpi_status_size) :: muster
5187 !      include 'COMMON.SETUP'
5188       logical :: tout,flag
5189       real(kind=8) :: tstart,tend1
5190       real(kind=8),parameter :: twait=600.0d0
5191       integer :: ihalt,man,ierr
5192
5193 !  find an available soldier
5194        tout=.false.
5195        flag=.false.
5196        tstart=MPI_WTIME()
5197        do while(.not. (flag .or. tout))
5198          call MPI_IPROBE(mpi_any_source,idint,CG_COMM,flag, &
5199                   muster,ierr)
5200          tend1=MPI_WTIME()
5201          if(tend1-tstart.gt.twait .and. ihalt.eq.1) tout=.true.
5202 !_error         if(tend1-tstart.gt.twait) tout=.true.
5203        enddo
5204        if (tout) then 
5205          write(iout,*) 'ERROR = timeout for recv ',tend1-tstart
5206          call flush(iout)
5207          return
5208        endif
5209        man=muster(mpi_source)        
5210
5211 !timeout         call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,
5212 !timeout     *                 CG_COMM,muster,ierr)
5213 !        print *, ' receiving output from start # ',ind(1)
5214 !t         print *,'receiving ',MPI_WTIME()
5215 !timeout         man=muster(mpi_source)
5216          call mpi_recv(ind,9,mpi_integer,man,idint,&
5217                        CG_COMM,muster,ierr)
5218 !timeout
5219 !  receive final energies and variables
5220          call mpi_recv(eout,1,mpi_double_precision,&
5221                        man,idreal,CG_COMM,muster,ierr)
5222 !         print *,eout 
5223 #ifdef CO_BIAS
5224          call mpi_recv(co,1,mpi_double_precision,&
5225                        man,idreal,CG_COMM,muster,ierr)
5226          write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co
5227 #endif
5228          call mpi_recv(xout,nvar,mpi_double_precision,&
5229                        man,idreal,CG_COMM,muster,ierr)
5230 !         print *,nvar , ierr
5231          if(vdisulf) nss=ind(6)
5232          if(vdisulf.and.nss.ne.0) then
5233           call mpi_recv(ihpb,nss,mpi_integer,&
5234                        man,idint,CG_COMM,muster,ierr)         
5235           call mpi_recv(jhpb,nss,mpi_integer,&
5236                        man,idint,CG_COMM,muster,ierr)
5237          endif
5238 !  halt soldier
5239        if(ihalt.eq.1) then 
5240 !        print *,'sending halt to ',man
5241         write(iout,*) 'sending halt to ',man
5242         info(1)=0
5243         call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,ierr)
5244        endif
5245       return
5246       end subroutine recv
5247 !-----------------------------------------------------------------------------
5248       subroutine history_append
5249
5250 !      implicit real*8 (a-h,o-z)
5251 !      include 'DIMENSIONS'
5252 !      include 'COMMON.IOUNITS'
5253
5254 #if defined(AIX) || defined(PGI)
5255        open(icsa_history,file=csa_history,position="append")
5256 #else
5257        open(icsa_history,file=csa_history,access="append")
5258 #endif
5259       return
5260       end subroutine history_append
5261 !-----------------------------------------------------------------------------
5262       subroutine alloc_CSA_arrays
5263
5264       use energy_data, only: ns
5265
5266       mxgr=2*nres
5267
5268       if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3))
5269 ! commom.bank
5270 !      common/varin/
5271 !el      allocate(dihang_in(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
5272       allocate(dihang_in(mxang,nres,mxch,5000)) !(mxang,maxres,mxch,mxio)
5273       allocate(nss_in(mxio)) !(mxio)
5274       allocate(iss_in(ns,mxio),jss_in(ns,mxio)) !(maxss,mxio)
5275 !      common/minvar/
5276       allocate(dihang(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
5277       allocate(rmsn(mxio),pncn(mxio)) !(mxio)
5278       allocate(etot(mxio)) !(mxio)
5279       allocate(nss_out(mxio)) !(mxio)
5280       allocate(iss_out(ns,mxio),jss_out(ns,mxio)) !(maxss,mxio)
5281 !      common/bank/
5282       allocate(rvar(mxang,nres,mxch,mxio),bvar(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio)
5283       allocate(bene(mxio),rene(mxio),brmsn(mxio),rrmsn(mxio))
5284       allocate(bpncn(mxio),rpncn(mxio)) !(mxio)
5285       allocate(ibank(mxio),is(mxio),jbank(mxio)) !(mxio)
5286       allocate(dij(mxio,mxio)) !(mxio,mxio)
5287 !      common/bank_disulfid/ 
5288       allocate(bvar_nss(mxio),bvar_ns(mxio)) !(mxio)
5289       allocate(bvar_s(ns,mxio)) !(maxss,mxio)
5290       allocate(bvar_ss(2,ns,mxio)) !(2,maxss,mxio) 
5291 !      common/mvstat/
5292       allocate(movenx(mxio),movernx(mxio)) !(mxio)
5293       allocate(nstatnx(0:mxmv,3),nstatnx_tot(0:mxmv,3)) !(0:mxmv,3)
5294       allocate(indb(mxio,9)) !(mxio,9)
5295       allocate(parent(3,mxio)) !(3,mxio)
5296 !      common/send2/
5297       allocate(isend2(mxio)) !(mxio)
5298       allocate(iff_in(nres,mxio2)) !(maxres,mxio2)
5299       allocate(dihang_in2(mxang,nres,mxch,mxio2)) !(mxang,maxres,mxch,mxio2)
5300       allocate(idata(5,mxio)) !(5,mxio)
5301 ! common.csa
5302 !      common/alphaa/
5303       allocate(ngroup(mxgr)) !(mxgr)
5304       allocate(igroup(3,mxang,mxgr)) !(3,mxang,mxgr)
5305 ! common.distfit
5306 !      COMMON /frag/
5307       allocate(bvar_frag(mxio,6)) !(mxio,6)
5308       allocate(hvar_frag(mxio,3),lvar_frag(mxio,3),svar_frag(mxio,3)) !(mxio,3)
5309       allocate(avar_frag(mxio,5)) !(mxio,5)
5310 ! commom.hairpin
5311 !      common /spinka/
5312       allocate(nharp_seed(nseed),nharp_use(nseed))      !(max_seed)
5313       allocate(iharp_seed(4,nres/3,nseed))      !(4,maxres/3,max_seed)
5314       allocate(iharp_use(0:4,nres/3,nseed))     !(0:4,maxres/3,max_seed)
5315
5316       return
5317       end subroutine alloc_CSA_arrays
5318 !-----------------------------------------------------------------------------
5319 !-----------------------------------------------------------------------------
5320       end module csa