wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / test.F
1       subroutine test
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4 #ifdef MPI
5       include 'mpif.h'
6 #endif
7       include 'COMMON.GEO'
8       include 'COMMON.VAR'
9       include 'COMMON.INTERACT'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.DISTFIT'
12       include 'COMMON.SBRIDGE'
13       include 'COMMON.CONTROL'
14       include 'COMMON.FFIELD'
15       include 'COMMON.MINIM'
16       include 'COMMON.CHAIN'
17       double precision time0,time1
18       double precision energy(0:n_ene),ee
19       double precision var(maxvar),var1(maxvar)
20       integer j1,j2
21       logical debug,accepted
22       debug=.true.
23       
24
25       call geom_to_var(nvar,var1)
26       call chainbuild
27       call etotal(energy(0))
28       etot=energy(0)
29       call rmsd(rms)
30       write(iout,*) 'etot=',0,etot,rms
31       call secondary2(.false.)
32
33       call write_pdb(0,'first structure',etot)
34
35       j1=13
36       j2=21
37       da=180.0*deg2rad
38
39
40
41        temp=3000.0d0
42        betbol=1.0D0/(1.9858D-3*temp)
43        jr=iran_num(j1,j2)
44        d=ran_number(-pi,pi)
45 c       phi(jr)=pinorm(phi(jr)+d)
46        call chainbuild
47        call etotal(energy(0))
48        etot0=energy(0)
49        call rmsd(rms)
50        write(iout,*) 'etot=',1,etot0,rms
51        call write_pdb(1,'perturb structure',etot0)
52
53       do i=2,500,2
54        jr=iran_num(j1,j2)
55        d=ran_number(-da,da)
56        phiold=phi(jr)
57        phi(jr)=pinorm(phi(jr)+d)
58        call chainbuild
59        call etotal(energy(0))
60        etot=energy(0)
61
62        if (etot.lt.etot0) then 
63           accepted=.true.
64        else
65           accepted=.false.
66           xxr=ran_number(0.0D0,1.0D0)
67           xxh=betbol*(etot-etot0)
68           if (xxh.lt.50.0D0) then
69             xxh=dexp(-xxh)
70             if (xxh.gt.xxr) accepted=.true. 
71           endif
72        endif
73        accepted=.true.
74 c       print *,etot0,etot,accepted
75        if (accepted) then 
76           etot0=etot
77           call rmsd(rms)
78           write(iout,*) 'etot=',i,etot,rms
79           call write_pdb(i,'MC structure',etot)
80 c minimize
81 c        call geom_to_var(nvar,var1)
82         call sc_move(2,nres-1,1,10d0,nft_sc,etot)
83         call geom_to_var(nvar,var)
84         call minimize(etot,var,iretcode,nfun)
85         write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
86         call var_to_geom(nvar,var)
87         call chainbuild
88         call rmsd(rms)
89         write(iout,*) 'etot mcm=',i,etot,rms
90         call write_pdb(i+1,'MCM structure',etot)
91         call var_to_geom(nvar,var1)
92 c --------
93        else
94           phi(jr)=phiold
95        endif
96       enddo
97
98 c minimize
99 c       call sc_move(2,nres-1,1,10d0,nft_sc,etot)
100 c       call geom_to_var(nvar,var)
101 c
102 c       call chainbuild        
103 c       call write_pdb(998 ,'sc min',etot)
104 c
105 c       call minimize(etot,var,iretcode,nfun)
106 c       write(iout,*)'------------------------------------------------'
107 c       write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun
108 c      
109 c       call var_to_geom(nvar,var)
110 c       call chainbuild        
111 c       call write_pdb(999,'full min',etot)
112
113
114       return
115       end
116
117
118       subroutine test_n16
119       implicit real*8 (a-h,o-z)
120       include 'DIMENSIONS'
121 #ifdef MPI
122       include 'mpif.h'
123 #endif
124       include 'COMMON.GEO'
125       include 'COMMON.VAR'
126       include 'COMMON.INTERACT'
127       include 'COMMON.IOUNITS'
128       include 'COMMON.DISTFIT'
129       include 'COMMON.SBRIDGE'
130       include 'COMMON.CONTROL'
131       include 'COMMON.FFIELD'
132       include 'COMMON.MINIM'
133       include 'COMMON.CHAIN'
134       double precision time0,time1
135       double precision energy(0:n_ene),ee
136       double precision var(maxvar),var1(maxvar)
137       integer jdata(5)
138       logical debug
139       debug=.true.
140
141 c
142       call geom_to_var(nvar,var1)
143       call chainbuild
144       call etotal(energy(0))
145       etot=energy(0)
146       write(iout,*) nnt,nct,etot
147       call write_pdb(1,'first structure',etot)
148       call secondary2(.true.)
149
150       do i=1,4
151         jdata(i)=bfrag(i,2)
152       enddo
153
154       DO ij=1,4
155        ieval=0
156        jdata(5)=ij
157        call var_to_geom(nvar,var1)
158        write(iout,*) 'N16 test',(jdata(i),i=1,5)
159        call beta_slide(jdata(1),jdata(2),jdata(3),jdata(4),jdata(5)
160      &               ,ieval,ij) 
161        call geom_to_var(nvar,var)       
162
163       if (minim) then
164 #ifdef MPI
165        time0=MPI_WTIME()
166 #else
167        time0=tcpu()
168 #endif
169        call minimize(etot,var,iretcode,nfun)
170        write(iout,*)'------------------------------------------------'
171        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
172      &  '+ DIST eval',ieval
173       
174 #ifdef MPI
175        time1=MPI_WTIME()
176 #else
177        time1=tcpu()
178 #endif
179        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
180      &         nfun/(time1-time0),' eval/s'
181
182        call var_to_geom(nvar,var)
183        call chainbuild        
184        call write_pdb(ij*100+99,'full min',etot)
185       endif
186
187
188       ENDDO
189
190       return
191       end
192
193
194       subroutine test_local
195       implicit real*8 (a-h,o-z)
196       include 'DIMENSIONS'
197       include 'COMMON.GEO'
198       include 'COMMON.VAR'
199       include 'COMMON.INTERACT'
200       include 'COMMON.IOUNITS'
201       double precision time0,time1
202       double precision energy(0:n_ene),ee
203       double precision varia(maxvar)
204 c
205       call chainbuild
206 c      call geom_to_var(nvar,varia)
207       call write_pdb(1,'first structure',0d0)
208
209       call etotal(energy(0))
210       etot=energy(0)
211       write(iout,*) nnt,nct,etot
212
213       write(iout,*) 'calling sc_move'
214       call sc_move(nnt,nct,5,10d0,nft_sc,etot)
215       write(iout,*) nft_sc,etot
216       call write_pdb(2,'second structure',etot)
217
218       write(iout,*) 'calling local_move'
219       call local_move_init(.false.)
220       call local_move(24,29,20d0,50d0)     
221       call chainbuild
222       call write_pdb(3,'third structure',etot)
223
224       write(iout,*) 'calling sc_move'
225       call sc_move(24,29,5,10d0,nft_sc,etot)
226       write(iout,*) nft_sc,etot
227       call write_pdb(2,'last structure',etot)
228
229
230       return
231       end
232
233       subroutine test_sc
234       implicit real*8 (a-h,o-z)
235       include 'DIMENSIONS'
236       include 'COMMON.GEO'
237       include 'COMMON.VAR'
238       include 'COMMON.INTERACT'
239       include 'COMMON.IOUNITS'
240       double precision time0,time1
241       double precision energy(0:n_ene),ee
242       double precision varia(maxvar)
243 c
244       call chainbuild
245 c      call geom_to_var(nvar,varia)
246       call write_pdb(1,'first structure',0d0)
247
248       call etotal(energy(0))
249       etot=energy(0)
250       write(iout,*) nnt,nct,etot
251
252       write(iout,*) 'calling sc_move'
253
254       call sc_move(nnt,nct,5,10d0,nft_sc,etot)
255       write(iout,*) nft_sc,etot
256       call write_pdb(2,'second structure',etot)
257
258       write(iout,*) 'calling sc_move 2nd time'
259
260       call sc_move(nnt,nct,5,1d0,nft_sc,etot)
261       write(iout,*) nft_sc,etot
262       call write_pdb(3,'last structure',etot)
263       return
264       end
265 c--------------------------------------------------------
266       subroutine bgrow(bstrand,nbstrand,in,ind,new)
267       implicit real*8 (a-h,o-z)
268       include 'DIMENSIONS'
269       include 'COMMON.CHAIN'
270       integer bstrand(maxres/3,6)
271
272       ishift=iabs(bstrand(in,ind+4)-new)
273
274       print *,'bgrow',bstrand(in,ind+4),new,ishift
275
276       bstrand(in,ind)=new
277
278       if(ind.eq.1)then
279         bstrand(nbstrand,5)=bstrand(nbstrand,1)
280         do i=1,nbstrand-1
281           IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
282           if (bstrand(i,5).lt.bstrand(i,6)) then 
283             bstrand(i,5)=bstrand(i,5)-ishift
284           else
285             bstrand(i,5)=bstrand(i,5)+ishift
286           endif
287           ENDIF
288         enddo
289       else
290         bstrand(nbstrand,6)=bstrand(nbstrand,2)
291         do i=1,nbstrand-1
292           IF (bstrand(nbstrand,3).eq.bstrand(i,3)) THEN
293           if (bstrand(i,6).lt.bstrand(i,5)) then 
294             bstrand(i,6)=bstrand(i,6)-ishift
295           else
296             bstrand(i,6)=bstrand(i,6)+ishift
297           endif
298           ENDIF
299         enddo
300       endif
301
302
303       return
304       end
305
306
307 c------------------------------------------
308       subroutine test11
309       implicit real*8 (a-h,o-z)
310       include 'DIMENSIONS'
311 #ifdef MPI
312       include 'mpif.h'
313 #endif
314       include 'COMMON.GEO'
315       include 'COMMON.CHAIN'
316       include 'COMMON.IOUNITS'
317       include 'COMMON.VAR'
318       include 'COMMON.CONTROL'
319       include 'COMMON.SBRIDGE'
320       include 'COMMON.FFIELD'
321       include 'COMMON.MINIM'
322 c
323       include 'COMMON.DISTFIT'       
324       integer if(20,maxres),nif,ifa(20)
325       integer ibc(0:maxres,0:maxres),istrand(20)
326       integer ibd(maxres),ifb(10,2),nifb,lifb(10),lifb0
327       integer itmp(20,maxres)
328       double precision time0,time1
329       double precision energy(0:n_ene),ee
330       double precision varia(maxvar),vorg(maxvar)
331 c
332       logical debug,ltest,usedbfrag(maxres/3)
333       character*50 linia
334 c
335       integer betasheet(maxres),ibetasheet(maxres),nbetasheet
336       integer bstrand(maxres/3,6),nbstrand
337
338 c------------------------ 
339
340       debug=.true.
341 c------------------------
342       nbstrand=0
343       nbetasheet=0
344       do i=1,nres
345         betasheet(i)=0
346         ibetasheet(i)=0
347       enddo
348       call geom_to_var(nvar,vorg)
349       call secondary2(debug)
350
351       if (nbfrag.le.1) return
352
353       do i=1,nbfrag
354          usedbfrag(i)=.false.
355       enddo
356
357
358       nbetasheet=nbetasheet+1
359       nbstrand=2
360       bstrand(1,1)=bfrag(1,1)
361       bstrand(1,2)=bfrag(2,1)
362       bstrand(1,3)=nbetasheet
363       bstrand(1,4)=1
364       bstrand(1,5)=bfrag(1,1)
365       bstrand(1,6)=bfrag(2,1)
366       do i=bfrag(1,1),bfrag(2,1)
367         betasheet(i)=nbetasheet
368         ibetasheet(i)=1
369       enddo
370 c
371       bstrand(2,1)=bfrag(3,1)
372       bstrand(2,2)=bfrag(4,1)
373       bstrand(2,3)=nbetasheet
374       bstrand(2,5)=bfrag(3,1)
375       bstrand(2,6)=bfrag(4,1)
376
377       if (bfrag(3,1).le.bfrag(4,1)) then
378         bstrand(2,4)=2
379         do i=bfrag(3,1),bfrag(4,1)
380           betasheet(i)=nbetasheet
381           ibetasheet(i)=2
382         enddo
383       else
384         bstrand(2,4)=-2
385         do i=bfrag(4,1),bfrag(3,1)
386           betasheet(i)=nbetasheet
387           ibetasheet(i)=2
388         enddo
389       endif
390
391       iused_nbfrag=1
392
393       do while (iused_nbfrag.ne.nbfrag)
394
395       do j=2,nbfrag
396        
397         IF (.not.usedbfrag(j)) THEN
398
399         write (*,*) j,(bfrag(i,j),i=1,4)
400         do jk=6,1,-1
401          write (*,'(i4,a3,10i4)') jk,'B',(bstrand(i,jk),i=1,nbstrand)
402         enddo
403         write (*,*) '------------------'
404
405
406         if (bfrag(3,j).le.bfrag(4,j)) then 
407          do i=bfrag(3,j),bfrag(4,j)
408           if(betasheet(i).eq.nbetasheet) then
409             in=ibetasheet(i)
410             do k=bfrag(3,j),bfrag(4,j)
411               betasheet(k)=nbetasheet
412               ibetasheet(k)=in
413             enddo
414             nbstrand=nbstrand+1
415             usedbfrag(j)=.true.
416             iused_nbfrag=iused_nbfrag+1
417             do k=bfrag(1,j),bfrag(2,j)
418               betasheet(k)=nbetasheet
419               ibetasheet(k)=nbstrand
420             enddo
421             if (bstrand(in,4).lt.0) then 
422               bstrand(nbstrand,1)=bfrag(2,j)
423               bstrand(nbstrand,2)=bfrag(1,j)
424               bstrand(nbstrand,3)=nbetasheet
425               bstrand(nbstrand,4)=-nbstrand
426               bstrand(nbstrand,5)=bstrand(nbstrand,1)
427               bstrand(nbstrand,6)=bstrand(nbstrand,2)
428               if(bstrand(in,1).lt.bfrag(4,j)) then
429                  call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
430               else
431                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+
432      &               (bstrand(in,5)-bfrag(4,j))
433               endif
434               if(bstrand(in,2).gt.bfrag(3,j)) then
435                  call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
436               else
437                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-
438      &               (-bstrand(in,6)+bfrag(3,j))                 
439               endif
440             else
441               bstrand(nbstrand,1)=bfrag(1,j)
442               bstrand(nbstrand,2)=bfrag(2,j)
443               bstrand(nbstrand,3)=nbetasheet
444               bstrand(nbstrand,4)=nbstrand
445               bstrand(nbstrand,5)=bstrand(nbstrand,1)
446               bstrand(nbstrand,6)=bstrand(nbstrand,2)
447               if(bstrand(in,1).gt.bfrag(3,j)) then
448                  call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
449               else
450                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-
451      &                (-bstrand(in,5)+bfrag(3,j))
452               endif
453               if(bstrand(in,2).lt.bfrag(4,j)) then
454                  call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
455               else
456                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+
457      &              (bstrand(in,6)-bfrag(4,j))
458               endif
459             endif
460             goto 11
461           endif
462           if(betasheet(bfrag(1,j)+i-bfrag(3,j)).eq.nbetasheet) then
463             in=ibetasheet(bfrag(1,j)+i-bfrag(3,j))
464             do k=bfrag(1,j),bfrag(2,j)
465               betasheet(k)=nbetasheet
466               ibetasheet(k)=in
467             enddo
468             nbstrand=nbstrand+1
469             usedbfrag(j)=.true.
470             iused_nbfrag=iused_nbfrag+1
471             do k=bfrag(3,1),bfrag(4,1)
472               betasheet(k)=nbetasheet
473               ibetasheet(k)=nbstrand
474             enddo
475             if (bstrand(in,4).lt.0) then 
476               bstrand(nbstrand,1)=bfrag(4,j)
477               bstrand(nbstrand,2)=bfrag(3,j)
478               bstrand(nbstrand,3)=nbetasheet
479               bstrand(nbstrand,4)=-nbstrand
480               bstrand(nbstrand,5)=bstrand(nbstrand,1)
481               bstrand(nbstrand,6)=bstrand(nbstrand,2)
482               if(bstrand(in,1).lt.bfrag(2,j)) then
483                  call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
484               else
485                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+
486      &              (bstrand(in,5)-bfrag(2,j))
487               endif
488               if(bstrand(in,2).gt.bfrag(1,j)) then
489                  call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
490               else
491                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-
492      &              (-bstrand(in,6)+bfrag(1,j))
493               endif
494             else
495               bstrand(nbstrand,1)=bfrag(3,j)
496               bstrand(nbstrand,2)=bfrag(4,j)
497               bstrand(nbstrand,3)=nbetasheet
498               bstrand(nbstrand,4)=nbstrand
499               bstrand(nbstrand,5)=bstrand(nbstrand,1)
500               bstrand(nbstrand,6)=bstrand(nbstrand,2)
501               if(bstrand(in,1).gt.bfrag(1,j)) then
502                  call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
503               else
504                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-
505      &              (-bstrand(in,5)+bfrag(1,j))
506               endif
507               if(bstrand(in,2).lt.bfrag(2,j)) then
508                  call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
509               else
510                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+
511      &               (bstrand(in,6)-bfrag(2,j))
512               endif
513             endif
514             goto 11
515           endif
516          enddo
517         else
518          do i=bfrag(4,j),bfrag(3,j)
519           if(betasheet(i).eq.nbetasheet) then
520             in=ibetasheet(i)
521             do k=bfrag(4,j),bfrag(3,j)
522               betasheet(k)=nbetasheet
523               ibetasheet(k)=in
524             enddo
525             nbstrand=nbstrand+1
526             usedbfrag(j)=.true.
527             iused_nbfrag=iused_nbfrag+1
528             do k=bfrag(1,j),bfrag(2,j)
529               betasheet(k)=nbetasheet
530               ibetasheet(k)=nbstrand
531             enddo
532             if (bstrand(in,4).lt.0) then 
533               bstrand(nbstrand,1)=bfrag(1,j)
534               bstrand(nbstrand,2)=bfrag(2,j)
535               bstrand(nbstrand,3)=nbetasheet
536               bstrand(nbstrand,4)=nbstrand
537               bstrand(nbstrand,5)=bstrand(nbstrand,1)
538               bstrand(nbstrand,6)=bstrand(nbstrand,2)
539               if(bstrand(in,1).lt.bfrag(3,j)) then
540                  call bgrow(bstrand,nbstrand,in,1,bfrag(3,j))
541               else
542                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-
543      &              (bstrand(in,5)-bfrag(3,j))
544               endif
545               if(bstrand(in,2).gt.bfrag(4,j)) then
546                  call bgrow(bstrand,nbstrand,in,2,bfrag(4,j))
547               else
548                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+
549      &              (-bstrand(in,6)+bfrag(4,j))
550               endif
551             else
552               bstrand(nbstrand,1)=bfrag(2,j)
553               bstrand(nbstrand,2)=bfrag(1,j)
554               bstrand(nbstrand,3)=nbetasheet
555               bstrand(nbstrand,4)=-nbstrand
556               bstrand(nbstrand,5)=bstrand(nbstrand,1)
557               bstrand(nbstrand,6)=bstrand(nbstrand,2)
558               if(bstrand(in,1).gt.bfrag(4,j)) then
559                  call bgrow(bstrand,nbstrand,in,1,bfrag(4,j))
560               else
561                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+
562      &              (-bstrand(in,5)+bfrag(4,j))
563               endif
564               if(bstrand(in,2).lt.bfrag(3,j)) then
565                  call bgrow(bstrand,nbstrand,in,2,bfrag(3,j))
566               else
567                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-
568      &             (bstrand(in,6)-bfrag(3,j))
569               endif
570             endif
571             goto 11
572           endif
573           if(betasheet(bfrag(2,j)-i+bfrag(4,j)).eq.nbetasheet) then
574             in=ibetasheet(bfrag(2,j)-i+bfrag(4,j))
575             do k=bfrag(1,j),bfrag(2,j)
576               betasheet(k)=nbetasheet
577               ibetasheet(k)=in
578             enddo
579             nbstrand=nbstrand+1
580             usedbfrag(j)=.true.
581             iused_nbfrag=iused_nbfrag+1
582             do k=bfrag(4,j),bfrag(3,j)  
583               betasheet(k)=nbetasheet
584               ibetasheet(k)=nbstrand
585             enddo
586             if (bstrand(in,4).lt.0) then 
587               bstrand(nbstrand,1)=bfrag(4,j)
588               bstrand(nbstrand,2)=bfrag(3,j)
589               bstrand(nbstrand,3)=nbetasheet
590               bstrand(nbstrand,4)=nbstrand
591               bstrand(nbstrand,5)=bstrand(nbstrand,1)
592               bstrand(nbstrand,6)=bstrand(nbstrand,2)
593               if(bstrand(in,1).lt.bfrag(2,j)) then
594                  call bgrow(bstrand,nbstrand,in,1,bfrag(2,j))
595               else
596                  bstrand(nbstrand,5)=bstrand(nbstrand,5)-
597      &             (bstrand(in,5)-bfrag(2,j))
598               endif
599               if(bstrand(in,2).gt.bfrag(1,j)) then 
600                  call bgrow(bstrand,nbstrand,in,2,bfrag(1,j))
601               else
602                  bstrand(nbstrand,6)=bstrand(nbstrand,6)+
603      &             (-bstrand(in,6)+bfrag(1,j))
604               endif
605             else
606               bstrand(nbstrand,1)=bfrag(3,j)
607               bstrand(nbstrand,2)=bfrag(4,j)
608               bstrand(nbstrand,3)=nbetasheet
609               bstrand(nbstrand,4)=-nbstrand
610               bstrand(nbstrand,5)=bstrand(nbstrand,1)
611               bstrand(nbstrand,6)=bstrand(nbstrand,2)
612               if(bstrand(in,1).gt.bfrag(1,j)) then
613                  call bgrow(bstrand,nbstrand,in,1,bfrag(1,j))
614               else
615                  bstrand(nbstrand,5)=bstrand(nbstrand,5)+
616      &               (-bstrand(in,5)+bfrag(1,j))
617               endif
618               if(bstrand(in,2).lt.bfrag(2,j)) then
619                  call bgrow(bstrand,nbstrand,in,2,bfrag(2,j))
620               else
621                  bstrand(nbstrand,6)=bstrand(nbstrand,6)-
622      &                (bstrand(in,6)-bfrag(2,j))
623               endif
624             endif
625             goto 11
626           endif
627          enddo
628         endif
629
630
631
632         ENDIF
633       enddo
634
635       j=2
636       do while (usedbfrag(j))      
637         j=j+1      
638       enddo
639    
640       nbstrand=nbstrand+1
641       nbetasheet=nbetasheet+1
642       bstrand(nbstrand,1)=bfrag(1,j)
643       bstrand(nbstrand,2)=bfrag(2,j)
644       bstrand(nbstrand,3)=nbetasheet
645       bstrand(nbstrand,5)=bfrag(1,j)
646       bstrand(nbstrand,6)=bfrag(2,j)
647       
648       bstrand(nbstrand,4)=nbstrand
649       do i=bfrag(1,j),bfrag(2,j)
650           betasheet(i)=nbetasheet
651           ibetasheet(i)=nbstrand
652       enddo
653 c
654       nbstrand=nbstrand+1
655       bstrand(nbstrand,1)=bfrag(3,j)
656       bstrand(nbstrand,2)=bfrag(4,j)
657       bstrand(nbstrand,3)=nbetasheet
658       bstrand(nbstrand,5)=bfrag(3,j)
659       bstrand(nbstrand,6)=bfrag(4,j)
660
661       if (bfrag(3,j).le.bfrag(4,j)) then
662         bstrand(nbstrand,4)=nbstrand
663         do i=bfrag(3,j),bfrag(4,j)
664           betasheet(i)=nbetasheet
665           ibetasheet(i)=nbstrand
666         enddo
667       else
668         bstrand(nbstrand,4)=-nbstrand
669         do i=bfrag(4,j),bfrag(3,j)
670           betasheet(i)=nbetasheet
671           ibetasheet(i)=nbstrand
672         enddo
673       endif
674
675       iused_nbfrag=iused_nbfrag+1
676       usedbfrag(j)=.true.
677
678
679   11  continue
680         do jk=6,1,-1
681          write (*,'(i4,a3,10i4)') jk,'A',(bstrand(i,jk),i=1,nbstrand)
682         enddo
683
684
685       enddo
686
687       do i=1,nres
688        if (betasheet(i).ne.0) write(*,*) i,betasheet(i),ibetasheet(i)
689       enddo
690       write(*,*)
691       do j=6,1,-1
692         write (*,'(i4,a3,10i4)') j,':',(bstrand(i,j),i=1,nbstrand)
693       enddo
694
695 c------------------------
696       nifb=0
697       do i=1,nbstrand
698         do j=i+1,nbstrand
699            if(iabs(bstrand(i,5)-bstrand(j,5)).le.5 .or.
700      &          iabs(bstrand(i,6)-bstrand(j,6)).le.5 ) then
701              nifb=nifb+1
702              ifb(nifb,1)=bstrand(i,4)
703              ifb(nifb,2)=bstrand(j,4)
704            endif
705         enddo
706       enddo
707
708       write(*,*)
709       do i=1,nifb
710          write (*,'(a3,20i4)') "ifb",i,ifb(i,1),ifb(i,2)
711       enddo
712
713       do i=1,nbstrand
714            ifa(i)=bstrand(i,4)
715       enddo
716       write (*,'(a3,20i4)') "ifa",(ifa(i),i=1,nbstrand)
717       
718       nif=iabs(bstrand(1,6)-bstrand(1,5))+1
719       do j=2,nbstrand
720        if (iabs(bstrand(j,6)-bstrand(j,5))+1.gt.nif)
721      &    nif=iabs(bstrand(j,6)-bstrand(j,5))+1
722       enddo
723      
724       write(*,*) nif
725       do i=1,nif
726         do j=1,nbstrand
727           if(j,i)=bstrand(j,6)+(i-1)*sign(1,bstrand(j,5)-bstrand(j,6))
728           if (if(j,i).gt.0) then
729             if(betasheet(if(j,i)).eq.0 .or. 
730      &          ibetasheet(if(j,i)).ne.iabs(bstrand(j,4))) if(j,i)=0  
731           else
732             if(j,i)=0
733           endif 
734         enddo
735         write(*,'(a3,10i4)') 'if ',(if(j,i),j=1,nbstrand)
736       enddo
737
738 c      read (inp,*) (ifa(i),i=1,4)    
739 c      do i=1,nres
740 c       read (inp,*,err=20,end=20) (if(j,i),j=1,4)
741 c      enddo
742 c 20   nif=i-1
743        stop
744 c------------------------
745
746       isa=4
747       is=2*isa-1
748       iconf=0
749 cccccccccccccccccccccccccccccccccc
750       DO ig=1,is**isa-1
751 cccccccccccccccccccccccccccccccccc
752
753          ii=ig
754          do j=1,is
755            istrand(is-j+1)=int(ii/is**(is-j))
756            ii=ii-istrand(is-j+1)*is**(is-j)
757          enddo  
758          ltest=.true.
759          do k=1,isa
760            istrand(k)=istrand(k)+1
761            if(istrand(k).gt.isa) istrand(k)=istrand(k)-2*isa-1
762          enddo
763          do k=1,isa
764            do l=1,isa
765             if(istrand(k).eq.istrand(l).and.k.ne.l.or.
766      &          istrand(k).eq.-istrand(l).and.k.ne.l) ltest=.false.
767            enddo
768          enddo
769
770          lifb0=1
771          do m=1,nifb
772            lifb(m)=0
773            do k=1,isa-1
774             if(
775      &         ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or.   
776      &         ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or.
777      &       -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or.
778      &       -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1))
779      &         lifb(m)=1
780            enddo
781          lifb0=lifb0*lifb(m)
782          enddo
783
784          if (mod(isa,2).eq.0) then
785           do k=isa/2+1,isa
786            if (istrand(k).eq.1) ltest=.false.
787           enddo
788          else
789           do k=(isa+1)/2+1,isa
790            if (istrand(k).eq.1) ltest=.false.
791           enddo          
792          endif
793
794          IF (ltest.and.lifb0.eq.1) THEN
795              iconf=iconf+1
796
797              call var_to_geom(nvar,vorg)
798
799              write (*,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
800              write (iout,'(i5,i10,10i3)') iconf,ig,(istrand(k),k=1,isa)
801              write (linia,'(10i3)') (istrand(k),k=1,isa)
802              
803       do i=1,nres
804         do j=1,nres
805          ibc(i,j)=0
806         enddo
807       enddo
808       
809
810       do i=1,4
811        if ( sign(1,istrand(i)).eq.sign(1,ifa(iabs(istrand(i)))) ) then
812         do j=1,nif
813          itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),j)
814         enddo         
815        else
816         do j=1,nif
817         itmp(iabs(istrand(i)),j)=if(iabs(ifa(iabs(istrand(i)))),nif-j+1)
818         enddo
819        endif
820       enddo
821      
822       do i=1,nif
823         write(*,*) (itmp(j,i),j=1,4)
824       enddo
825
826       do i=1,nif
827 c       ifa(1),ifa(2),ifa(3),ifa(4)
828 c       if(1,i),if(2,i),if(3,i),if(4,i)
829         do k=1,isa-1
830          ltest=.false.
831          do m=1,nifb
832            if(
833      &         ifb(m,1).eq.istrand(k).and.ifb(m,2).eq.istrand(k+1).or.   
834      &         ifb(m,2).eq.istrand(k).and.ifb(m,1).eq.istrand(k+1).or.
835      &       -ifb(m,1).eq.istrand(k).and.-ifb(m,2).eq.istrand(k+1).or.
836      &       -ifb(m,2).eq.istrand(k).and.-ifb(m,1).eq.istrand(k+1))
837      &   then 
838            ltest=.true.
839            goto 110
840          endif  
841          enddo
842   110     continue
843          if (ltest) then
844           ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-1
845          else
846           ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+1)),i))=-2
847          endif
848 c
849         if (k.lt.3)
850      &   ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+2)),i))=-3
851         if (k.lt.2)
852      &   ibc(itmp(iabs(istrand(k)),i),itmp(iabs(istrand(k+3)),i))=-4
853         enddo
854       enddo
855 c------------------------
856
857 c
858 c  freeze sec.elements 
859 c
860        do i=1,nres
861          mask(i)=1
862          mask_phi(i)=1
863          mask_theta(i)=1
864          mask_side(i)=1
865        enddo
866
867        do j=1,nbfrag
868         do i=bfrag(1,j),bfrag(2,j)
869          mask(i)=0
870          mask_phi(i)=0
871          mask_theta(i)=0
872         enddo
873         if (bfrag(3,j).le.bfrag(4,j)) then 
874          do i=bfrag(3,j),bfrag(4,j)
875           mask(i)=0
876           mask_phi(i)=0
877           mask_theta(i)=0
878          enddo
879         else
880          do i=bfrag(4,j),bfrag(3,j)
881           mask(i)=0
882           mask_phi(i)=0
883           mask_theta(i)=0
884          enddo
885         endif
886        enddo
887        do j=1,nhfrag
888         do i=hfrag(1,j),hfrag(2,j)
889          mask(i)=0
890          mask_phi(i)=0
891          mask_theta(i)=0
892         enddo
893        enddo
894        mask_r=.true.
895
896 c------------------------
897 c      generate constrains 
898 c
899        nhpb0=nhpb
900        call chainbuild                                                           
901        ind=0                                                                     
902        do i=1,nres-3                                                             
903          do j=i+3,nres                                                           
904           ind=ind+1                                                              
905           if ( ibc(i,j).eq.-1 .or. ibc(j,i).eq.-1) then                                           
906             d0(ind)=DIST(i,j)                                                     
907             w(ind)=10.0                                                           
908             nhpb=nhpb+1                                                           
909             ihpb(nhpb)=i                                                          
910             jhpb(nhpb)=j                                                          
911             forcon(nhpb)=10.0                                                     
912             dhpb(nhpb)=d0(ind)         
913           else if ( ibc(i,j).eq.-2 .or. ibc(j,i).eq.-2) then
914             d0(ind)=5.0
915             w(ind)=10.0                                                           
916             nhpb=nhpb+1                                                           
917             ihpb(nhpb)=i                                                          
918             jhpb(nhpb)=j
919             forcon(nhpb)=10.0                                                     
920             dhpb(nhpb)=d0(ind)                                                    
921           else if ( ibc(i,j).eq.-3 .or. ibc(j,i).eq.-3) then
922             d0(ind)=11.0
923             w(ind)=10.0                                                           
924             nhpb=nhpb+1                                                           
925             ihpb(nhpb)=i                                                          
926             jhpb(nhpb)=j
927             forcon(nhpb)=10.0                                                     
928             dhpb(nhpb)=d0(ind)                                                    
929           else if ( ibc(i,j).eq.-4 .or. ibc(j,i).eq.-4) then
930             d0(ind)=16.0
931             w(ind)=10.0                                                           
932             nhpb=nhpb+1                                                           
933             ihpb(nhpb)=i                                                          
934             jhpb(nhpb)=j
935             forcon(nhpb)=10.0                                                     
936             dhpb(nhpb)=d0(ind)                                                    
937           else if ( ibc(i,j).gt.0 ) then
938             d0(ind)=DIST(i,ibc(i,j))             
939             w(ind)=10.0                                                           
940             nhpb=nhpb+1                                                           
941             ihpb(nhpb)=i                                                          
942             jhpb(nhpb)=j
943             forcon(nhpb)=10.0                                                     
944             dhpb(nhpb)=d0(ind)         
945           else if ( ibc(j,i).gt.0 ) then
946             d0(ind)=DIST(ibc(j,i),j)             
947             w(ind)=10.0                                                           
948             nhpb=nhpb+1                                                           
949             ihpb(nhpb)=i                                                          
950             jhpb(nhpb)=j
951             forcon(nhpb)=10.0                                                     
952             dhpb(nhpb)=d0(ind)         
953           else
954             w(ind)=0.0
955           endif                                                                  
956           dd(ind)=d0(ind)
957          enddo                                                                   
958        enddo                                    
959        call hpb_partition
960 cd--------------------------
961
962       write(iout,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),
963      &        ibc(jhpb(i),ihpb(i)),' --',
964      &        ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
965
966 cd      nhpb=0
967 cd      goto 901
968 c
969 c
970       call contact_cp_min(varia,ifun,iconf,linia,debug)
971       if (minim) then
972 #ifdef MPI
973        time0=MPI_WTIME()
974 #else
975        time0=tcpu()
976 #endif
977        call minimize(etot,varia,iretcode,nfun)
978        write(iout,*)'------------------------------------------------'
979        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
980      &  '+ DIST eval',ifun
981       
982 #ifdef MPI
983        time1=MPI_WTIME()
984 #else
985        time1=tcpu()
986 #endif
987        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
988      &         nfun/(time1-time0),' eval/s'
989
990        write (linia,'(a10,10i3)') 'full_min',(istrand(k),k=1,isa)
991        call var_to_geom(nvar,varia)
992        call chainbuild 
993        call write_pdb(900+iconf,linia,etot)
994       endif
995        
996       call etotal(energy(0))
997       etot=energy(0)
998       call enerprint(energy(0))
999 cd      call intout      
1000 cd      call briefout(0,etot)
1001 cd      call secondary2(.true.)
1002
1003  901  CONTINUE 
1004 ctest      return
1005 cccccccccccccccccccccccccccccccccccc
1006       ENDIF
1007       ENDDO
1008 cccccccccccccccccccccccccccccccccccc
1009
1010       return
1011   10  write (iout,'(a)') 'Error reading test structure.'
1012       return
1013       end
1014 c--------------------------------------------------------
1015
1016       subroutine test3
1017       implicit real*8 (a-h,o-z)
1018       include 'DIMENSIONS'
1019 #ifdef MPI
1020       include 'mpif.h'
1021 #endif
1022       include 'COMMON.GEO'
1023       include 'COMMON.CHAIN'
1024       include 'COMMON.IOUNITS'
1025       include 'COMMON.VAR'
1026       include 'COMMON.CONTROL'
1027       include 'COMMON.SBRIDGE'
1028       include 'COMMON.FFIELD'
1029       include 'COMMON.MINIM'
1030 c
1031       include 'COMMON.DISTFIT'       
1032       integer if(3,maxres),nif
1033       integer ibc(maxres,maxres),istrand(20)
1034       integer ibd(maxres),ifb(10,2),nifb,lifb(10),lifb0
1035       double precision time0,time1
1036       double precision energy(0:n_ene),ee
1037       double precision varia(maxvar)
1038 c
1039       logical debug,ltest
1040       character*50 linia
1041 c
1042       do i=1,nres
1043        read (inp,*,err=20,end=20) if(1,i),if(2,i),if(3,i)
1044       enddo
1045  20   nif=i-1
1046       write (*,'(a4,3i5)') ('if =',if(1,i),if(2,i),if(3,i),
1047      &                                         i=1,nif)
1048
1049         
1050 c------------------------
1051       call secondary2(debug)
1052 c------------------------
1053       do i=1,nres
1054         do j=1,nres
1055          ibc(i,j)=0
1056         enddo
1057       enddo
1058
1059 c
1060 c  freeze sec.elements and store indexes for beta constrains
1061 c
1062        do i=1,nres
1063          mask(i)=1
1064          mask_phi(i)=1
1065          mask_theta(i)=1
1066          mask_side(i)=1
1067        enddo
1068
1069        do j=1,nbfrag
1070         do i=bfrag(1,j),bfrag(2,j)
1071          mask(i)=0
1072          mask_phi(i)=0
1073          mask_theta(i)=0
1074         enddo
1075         if (bfrag(3,j).le.bfrag(4,j)) then 
1076          do i=bfrag(3,j),bfrag(4,j)
1077           mask(i)=0
1078           mask_phi(i)=0
1079           mask_theta(i)=0
1080           ibc(bfrag(1,j)+i-bfrag(3,j),i)=-1
1081          enddo
1082         else
1083          do i=bfrag(4,j),bfrag(3,j)
1084           mask(i)=0
1085           mask_phi(i)=0
1086           mask_theta(i)=0
1087           ibc(bfrag(2,j)-i+bfrag(4,j),i)=-1
1088          enddo
1089         endif
1090        enddo
1091        do j=1,nhfrag
1092         do i=hfrag(1,j),hfrag(2,j)
1093          mask(i)=0
1094          mask_phi(i)=0
1095          mask_theta(i)=0
1096         enddo
1097        enddo
1098        mask_r=.true.
1099
1100         
1101 c ---------------- test --------------
1102        do i=1,nif
1103          if (ibc(if(1,i),if(2,i)).eq.-1) then
1104              ibc(if(1,i),if(2,i))=if(3,i)
1105              ibc(if(1,i),if(3,i))=if(2,i)
1106          else if (ibc(if(2,i),if(1,i)).eq.-1) then
1107              ibc(if(2,i),if(1,i))=0
1108              ibc(if(1,i),if(2,i))=if(3,i)
1109              ibc(if(1,i),if(3,i))=if(2,i)
1110          else
1111              ibc(if(1,i),if(2,i))=if(3,i)
1112              ibc(if(1,i),if(3,i))=if(2,i)
1113          endif
1114        enddo
1115
1116        do i=1,nres
1117         do j=1,nres
1118          if (ibc(i,j).ne.0)  write(*,'(3i5)') i,j,ibc(i,j)
1119         enddo
1120        enddo
1121 c------------------------
1122        call chainbuild                                                           
1123        ind=0                                                                     
1124        do i=1,nres-3                                                             
1125          do j=i+3,nres                                                           
1126           ind=ind+1                                                              
1127           if ( ibc(i,j).eq.-1 ) then                                           
1128             d0(ind)=DIST(i,j)                                                     
1129             w(ind)=10.0                                                           
1130             nhpb=nhpb+1                                                           
1131             ihpb(nhpb)=i                                                          
1132             jhpb(nhpb)=j                                                          
1133             forcon(nhpb)=10.0                                                     
1134             dhpb(nhpb)=d0(ind)                                                    
1135           else if ( ibc(i,j).gt.0 ) then
1136             d0(ind)=DIST(i,ibc(i,j))             
1137             w(ind)=10.0                                                           
1138             nhpb=nhpb+1                                                           
1139             ihpb(nhpb)=i                                                          
1140             jhpb(nhpb)=j
1141             forcon(nhpb)=10.0                                                     
1142             dhpb(nhpb)=d0(ind)         
1143           else if ( ibc(j,i).gt.0 ) then
1144             d0(ind)=DIST(ibc(j,i),j)             
1145             w(ind)=10.0                                                           
1146             nhpb=nhpb+1                                                           
1147             ihpb(nhpb)=i                                                          
1148             jhpb(nhpb)=j
1149             forcon(nhpb)=10.0                                                     
1150             dhpb(nhpb)=d0(ind)         
1151           else
1152             w(ind)=0.0
1153           endif                                                                  
1154          enddo                                                                   
1155        enddo                                    
1156        call hpb_partition
1157
1158 cd--------------------------
1159       write(*,'(i3,2i4,a3,2i4,f7.2)') (i,ibc(ihpb(i),jhpb(i)),
1160      &        ibc(jhpb(i),ihpb(i)),' --',
1161      &        ihpb(i),jhpb(i),dhpb(i),i=1,nhpb)
1162       
1163
1164       linia='dist'
1165       debug=.true.
1166       in_pdb=7
1167 c
1168       call contact_cp_min(varia,ieval,in_pdb,linia,debug)
1169       if (minim) then
1170 #ifdef MPI
1171        time0=MPI_WTIME()
1172 #else
1173        time0=tcpu()
1174 #endif
1175        call minimize(etot,varia,iretcode,nfun)
1176        write(iout,*)'------------------------------------------------'
1177        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
1178      &  '+ DIST eval',ieval
1179       
1180 #ifdef MPI
1181        time1=MPI_WTIME()
1182 #else
1183        time1=tcpu()
1184 #endif
1185        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
1186      &         nfun/(time1-time0),' eval/s'
1187
1188
1189        call var_to_geom(nvar,varia)
1190        call chainbuild        
1191        call write_pdb(999,'full min',etot)
1192       endif
1193        
1194       call etotal(energy(0))
1195       etot=energy(0)
1196       call enerprint(energy(0))
1197       call intout      
1198       call briefout(0,etot)
1199       call secondary2(.true.)
1200
1201       return
1202   10  write (iout,'(a)') 'Error reading test structure.'
1203       return
1204       end
1205
1206
1207
1208
1209       subroutine test__
1210       implicit real*8 (a-h,o-z)
1211       include 'DIMENSIONS'
1212 #ifdef MPI
1213       include 'mpif.h'
1214 #endif
1215       include 'COMMON.GEO'
1216       include 'COMMON.CHAIN'
1217       include 'COMMON.IOUNITS'
1218       include 'COMMON.VAR'
1219       include 'COMMON.CONTROL'
1220       include 'COMMON.SBRIDGE'
1221       include 'COMMON.FFIELD'
1222       include 'COMMON.MINIM'
1223 c
1224       include 'COMMON.DISTFIT'       
1225       integer if(2,2),ind
1226       integer iff(maxres)
1227       double precision time0,time1
1228       double precision energy(0:n_ene),ee
1229       double precision theta2(maxres),phi2(maxres),alph2(maxres),
1230      &              omeg2(maxres),
1231      &              theta1(maxres),phi1(maxres),alph1(maxres),
1232      &              omeg1(maxres)
1233       double precision varia(maxvar),varia2(maxvar)
1234 c
1235
1236
1237       read (inp,*,err=10,end=10) if(1,1),if(1,2),if(2,1),if(2,2)
1238       write (iout,'(a4,4i5)') 'if =',if(1,1),if(1,2),if(2,1),if(2,2)
1239       read (inp,*,err=10,end=10) (theta2(i),i=3,nres)                          
1240       read (inp,*,err=10,end=10) (phi2(i),i=4,nres)                            
1241       read (inp,*,err=10,end=10) (alph2(i),i=2,nres-1)                         
1242       read (inp,*,err=10,end=10) (omeg2(i),i=2,nres-1)                         
1243       do i=1,nres                                                               
1244         theta2(i)=deg2rad*theta2(i)                                               
1245         phi2(i)=deg2rad*phi2(i)                                                   
1246         alph2(i)=deg2rad*alph2(i)                                                 
1247         omeg2(i)=deg2rad*omeg2(i)                                                 
1248       enddo  
1249       do i=1,nres                                                               
1250         theta1(i)=theta(i)                                               
1251         phi1(i)=phi(i)                                                   
1252         alph1(i)=alph(i)                                                 
1253         omeg1(i)=omeg(i)                                                 
1254       enddo  
1255
1256       do i=1,nres
1257        mask(i)=1
1258       enddo
1259
1260       
1261 c------------------------
1262       do i=1,nres
1263          iff(i)=0
1264       enddo
1265       do j=1,2
1266         do i=if(j,1),if(j,2)
1267           iff(i)=1
1268         enddo
1269       enddo
1270
1271       call chainbuild
1272       call geom_to_var(nvar,varia)
1273       call write_pdb(1,'first structure',0d0)
1274
1275         call secondary(.true.)
1276         
1277         call secondary2(.true.)
1278
1279         do j=1,nbfrag
1280            if ( (bfrag(3,j).lt.bfrag(4,j) .or.
1281      &          bfrag(4,j)-bfrag(2,j).gt.4) .and. 
1282      &          bfrag(2,j)-bfrag(1,j).gt.3 ) then
1283              nn=nn+1
1284
1285              if (bfrag(3,j).lt.bfrag(4,j)) then
1286                write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') 
1287      &                     "select",bfrag(1,j)-1,"-",bfrag(2,j)-1
1288      &                         ,",",bfrag(3,j)-1,"-",bfrag(4,j)-1
1289              else
1290                write(iout,'(a6,i3,a1,i3,a1,i3,a1,i3)') 
1291      &                     "select",bfrag(1,j)-1,"-",bfrag(2,j)-1
1292      &                         ,",",bfrag(4,j)-1,"-",bfrag(3,j)-1
1293              endif
1294            endif
1295         enddo
1296
1297       do i=1,nres                                                               
1298         theta(i)=theta2(i)                                               
1299         phi(i)=phi2(i)                                                   
1300         alph(i)=alph2(i)                                                 
1301         omeg(i)=omeg2(i)                                                 
1302       enddo  
1303
1304       call chainbuild
1305       call geom_to_var(nvar,varia2)
1306       call write_pdb(2,'second structure',0d0)
1307
1308
1309
1310 c-------------------------------------------------------
1311
1312       ifun=-1
1313       call contact_cp(varia,varia2,iff,ifun,7)
1314       if (minim) then
1315 #ifdef MPI
1316        time0=MPI_WTIME()
1317 #else
1318        time0=tcpu()
1319 #endif
1320        call minimize(etot,varia,iretcode,nfun)
1321        write(iout,*)'------------------------------------------------'
1322        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
1323      &  '+ DIST eval',ifun
1324       
1325 #ifdef MPI
1326        time1=MPI_WTIME()
1327 #else
1328        time1=tcpu()
1329 #endif
1330        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
1331      &         nfun/(time1-time0),' eval/s'
1332
1333
1334        call var_to_geom(nvar,varia)
1335        call chainbuild        
1336        call write_pdb(999,'full min',etot)
1337       endif
1338        
1339       call etotal(energy(0))
1340       etot=energy(0)
1341       call enerprint(energy(0))
1342       call intout      
1343       call briefout(0,etot)
1344
1345       return
1346   10  write (iout,'(a)') 'Error reading test structure.'
1347       return
1348       end
1349
1350 c-------------------------------------------------
1351 c-------------------------------------------------
1352
1353       subroutine secondary(lprint)
1354       implicit real*8 (a-h,o-z)
1355       include 'DIMENSIONS'
1356       include 'COMMON.CHAIN'
1357       include 'COMMON.IOUNITS'
1358       include 'COMMON.DISTFIT'
1359
1360       integer ncont,icont(2,maxres*maxres/2),isec(maxres,3)
1361       logical lprint,not_done
1362       real dcont(maxres*maxres/2),d
1363       real rcomp /7.0/ 
1364       real rbeta /5.2/
1365       real ralfa /5.2/
1366       real r310 /6.6/
1367       double precision xpi(3),xpj(3)
1368
1369
1370
1371       call chainbuild
1372 cd      call write_pdb(99,'sec structure',0d0)
1373       ncont=0
1374       nbfrag=0
1375       nhfrag=0
1376       do i=1,nres
1377         isec(i,1)=0
1378         isec(i,2)=0
1379         isec(i,3)=0
1380       enddo
1381
1382       do i=2,nres-3
1383         do k=1,3        
1384           xpi(k)=0.5d0*(c(k,i-1)+c(k,i))
1385         enddo
1386         do j=i+2,nres
1387           do k=1,3
1388              xpj(k)=0.5d0*(c(k,j-1)+c(k,j))
1389           enddo
1390 cd        d = (c(1,i)-c(1,j))*(c(1,i)-c(1,j)) +
1391 cd     &         (c(2,i)-c(2,j))*(c(2,i)-c(2,j)) +
1392 cd     &         (c(3,i)-c(3,j))*(c(3,i)-c(3,j)) 
1393 cd          print *,'CA',i,j,d
1394           d =  (xpi(1)-xpj(1))*(xpi(1)-xpj(1)) +
1395      &         (xpi(2)-xpj(2))*(xpi(2)-xpj(2)) +
1396      &         (xpi(3)-xpj(3))*(xpi(3)-xpj(3)) 
1397           if ( d.lt.rcomp*rcomp) then
1398             ncont=ncont+1
1399             icont(1,ncont)=i
1400             icont(2,ncont)=j
1401             dcont(ncont)=sqrt(d)
1402           endif
1403         enddo
1404       enddo
1405       if (lprint) then
1406         write (iout,*)
1407         write (iout,'(a)') '#PP contact map distances:'
1408         do i=1,ncont
1409           write (iout,'(3i4,f10.5)') 
1410      &     i,icont(1,i),icont(2,i),dcont(i) 
1411         enddo
1412       endif
1413
1414 c finding parallel beta
1415 cd      write (iout,*) '------- looking for parallel beta -----------'
1416       nbeta=0
1417       nstrand=0
1418       do i=1,ncont
1419         i1=icont(1,i)
1420         j1=icont(2,i)
1421         if(dcont(i).le.rbeta .and. j1-i1.gt.4 .and.
1422      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
1423      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. 
1424      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. 
1425      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. 
1426      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
1427      &     ) then
1428           ii1=i1
1429           jj1=j1
1430 cd         write (iout,*) i1,j1,dcont(i)
1431           not_done=.true.
1432           do while (not_done)
1433            i1=i1+1
1434            j1=j1+1
1435             do j=1,ncont
1436               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)
1437      &              .and. dcont(j).le.rbeta .and.
1438      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
1439      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. 
1440      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. 
1441      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. 
1442      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
1443      &                            ) goto 5
1444             enddo
1445             not_done=.false.
1446   5         continue
1447 cd            write (iout,*) i1,j1,dcont(j),not_done
1448           enddo
1449           j1=j1-1
1450           i1=i1-1
1451           if (i1-ii1.gt.1) then
1452             ii1=max0(ii1-1,1)
1453             jj1=max0(jj1-1,1)
1454             nbeta=nbeta+1
1455             if(lprint)write(iout,*)'parallel beta',nbeta,ii1,i1,jj1,j1
1456
1457             nbfrag=nbfrag+1
1458             bfrag(1,nbfrag)=ii1
1459             bfrag(2,nbfrag)=i1
1460             bfrag(3,nbfrag)=jj1
1461             bfrag(4,nbfrag)=j1 
1462
1463             do ij=ii1,i1
1464              isec(ij,1)=isec(ij,1)+1
1465              isec(ij,1+isec(ij,1))=nbeta
1466             enddo
1467             do ij=jj1,j1
1468              isec(ij,1)=isec(ij,1)+1
1469              isec(ij,1+isec(ij,1))=nbeta
1470             enddo
1471
1472            if(lprint) then 
1473             nstrand=nstrand+1
1474             if (nbeta.le.9) then
1475               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
1476      &          "DefPropRes 'strand",nstrand,
1477      &          "' 'num = ",ii1-1,"..",i1-1,"'"
1478             else
1479               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
1480      &          "DefPropRes 'strand",nstrand,
1481      &          "' 'num = ",ii1-1,"..",i1-1,"'"
1482             endif
1483             nstrand=nstrand+1
1484             if (nbeta.le.9) then
1485               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
1486      &          "DefPropRes 'strand",nstrand,
1487      &          "' 'num = ",jj1-1,"..",j1-1,"'"
1488             else
1489               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
1490      &          "DefPropRes 'strand",nstrand,
1491      &          "' 'num = ",jj1-1,"..",j1-1,"'"
1492             endif
1493               write(12,'(a8,4i4)')
1494      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
1495            endif
1496           endif
1497         endif
1498       enddo
1499
1500 c finding antiparallel beta
1501 cd      write (iout,*) '--------- looking for antiparallel beta ---------'
1502
1503       do i=1,ncont
1504         i1=icont(1,i)
1505         j1=icont(2,i)
1506         if (dcont(i).le.rbeta.and.
1507      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
1508      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. 
1509      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. 
1510      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. 
1511      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
1512      &     ) then
1513           ii1=i1
1514           jj1=j1
1515 cd          write (iout,*) i1,j1,dcont(i)
1516
1517           not_done=.true.
1518           do while (not_done)
1519            i1=i1+1
1520            j1=j1-1
1521             do j=1,ncont
1522               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
1523      &      isec(i1,1).le.1.and.isec(j1,1).le.1.and.
1524      &    (isec(i1,2).ne.isec(j1,2).or.isec(i1,2)*isec(j1,2).eq.0).and. 
1525      &    (isec(i1,3).ne.isec(j1,3).or.isec(i1,3)*isec(j1,3).eq.0).and. 
1526      &    (isec(i1,2).ne.isec(j1,3).or.isec(i1,2)*isec(j1,3).eq.0).and. 
1527      &    (isec(i1,3).ne.isec(j1,2).or.isec(i1,3)*isec(j1,2).eq.0)
1528      &           .and. dcont(j).le.rbeta ) goto 6
1529             enddo
1530             not_done=.false.
1531   6         continue
1532 cd            write (iout,*) i1,j1,dcont(j),not_done
1533           enddo
1534           i1=i1-1
1535           j1=j1+1
1536           if (i1-ii1.gt.1) then
1537             if(lprint)write (iout,*)'antiparallel beta',
1538      &                   nbeta,ii1-1,i1,jj1,j1-1
1539
1540             nbfrag=nbfrag+1
1541             bfrag(1,nbfrag)=max0(ii1-1,1)
1542             bfrag(2,nbfrag)=i1
1543             bfrag(3,nbfrag)=jj1
1544             bfrag(4,nbfrag)=max0(j1-1,1) 
1545
1546             nbeta=nbeta+1
1547             iii1=max0(ii1-1,1)
1548             do ij=iii1,i1
1549              isec(ij,1)=isec(ij,1)+1
1550              isec(ij,1+isec(ij,1))=nbeta
1551             enddo
1552             jjj1=max0(j1-1,1)  
1553             do ij=jjj1,jj1
1554              isec(ij,1)=isec(ij,1)+1
1555              isec(ij,1+isec(ij,1))=nbeta
1556             enddo
1557
1558
1559            if (lprint) then
1560             nstrand=nstrand+1
1561             if (nstrand.le.9) then
1562               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
1563      &          "DefPropRes 'strand",nstrand,
1564      &          "' 'num = ",ii1-2,"..",i1-1,"'"
1565             else
1566               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
1567      &          "DefPropRes 'strand",nstrand,
1568      &          "' 'num = ",ii1-2,"..",i1-1,"'"
1569             endif
1570             nstrand=nstrand+1
1571             if (nstrand.le.9) then
1572               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
1573      &          "DefPropRes 'strand",nstrand,
1574      &          "' 'num = ",j1-2,"..",jj1-1,"'"
1575             else
1576               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
1577      &          "DefPropRes 'strand",nstrand,
1578      &          "' 'num = ",j1-2,"..",jj1-1,"'"
1579             endif
1580               write(12,'(a8,4i4)')
1581      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
1582            endif
1583           endif
1584         endif
1585       enddo
1586
1587       if (nstrand.gt.0.and.lprint) then
1588         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
1589         do i=2,nstrand
1590          if (i.le.9) then
1591           write(12,'(a9,i1,$)') " | strand",i
1592          else
1593           write(12,'(a9,i2,$)') " | strand",i
1594          endif
1595         enddo
1596         write(12,'(a1)') "'"
1597       endif
1598
1599        
1600 c finding alpha or 310 helix
1601
1602       nhelix=0
1603       do i=1,ncont
1604         i1=icont(1,i)
1605         j1=icont(2,i)
1606         if (j1.eq.i1+3.and.dcont(i).le.r310
1607      &     .or.j1.eq.i1+4.and.dcont(i).le.ralfa ) then
1608 cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,dcont(i)
1609 cd          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,dcont(i)
1610           ii1=i1
1611           jj1=j1
1612           if (isec(ii1,1).eq.0) then 
1613             not_done=.true.
1614           else
1615             not_done=.false.
1616           endif
1617           do while (not_done)
1618             i1=i1+1
1619             j1=j1+1
1620             do j=1,ncont
1621               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
1622             enddo
1623             not_done=.false.
1624   10        continue
1625 cd            write (iout,*) i1,j1,not_done
1626           enddo
1627           j1=j1-1
1628           if (j1-ii1.gt.4) then
1629             nhelix=nhelix+1
1630 cd            write (iout,*)'helix',nhelix,ii1,j1
1631
1632             nhfrag=nhfrag+1
1633             hfrag(1,nhfrag)=ii1
1634             hfrag(2,nhfrag)=max0(j1-1,1)
1635
1636             do ij=ii1,j1
1637              isec(ij,1)=-1
1638             enddo
1639            if (lprint) then
1640             write (iout,'(a6,i3,2i4)') "Helix",nhelix,ii1-1,j1-2
1641             if (nhelix.le.9) then
1642               write(12,'(a17,i1,a9,i3,a2,i3,a1)') 
1643      &          "DefPropRes 'helix",nhelix,
1644      &          "' 'num = ",ii1-1,"..",j1-2,"'"
1645             else
1646               write(12,'(a17,i2,a9,i3,a2,i3,a1)') 
1647      &          "DefPropRes 'helix",nhelix,
1648      &          "' 'num = ",ii1-1,"..",j1-2,"'"
1649             endif
1650            endif
1651           endif
1652         endif
1653       enddo
1654        
1655       if (nhelix.gt.0.and.lprint) then
1656         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
1657         do i=2,nhelix
1658          if (nhelix.le.9) then
1659           write(12,'(a8,i1,$)') " | helix",i
1660          else
1661           write(12,'(a8,i2,$)') " | helix",i
1662          endif
1663         enddo
1664         write(12,'(a1)') "'"
1665       endif
1666
1667       if (lprint) then
1668        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
1669        write(12,'(a20)') "XMacStand ribbon.mac"
1670       endif
1671
1672
1673       return
1674       end
1675 c----------------------------------------------------------------------------
1676
1677       subroutine write_pdb(npdb,titelloc,ee)
1678       implicit real*8 (a-h,o-z)
1679       include 'DIMENSIONS'
1680       include 'COMMON.IOUNITS'
1681       character*50 titelloc1                                                     
1682       character*(*) titelloc
1683       character*3 zahl   
1684       character*5 liczba5
1685       double precision ee
1686       integer npdb,ilen
1687       external ilen
1688
1689       titelloc1=titelloc
1690       lenpre=ilen(prefix)
1691       if (npdb.lt.1000) then
1692        call numstr(npdb,zahl)
1693        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
1694       else
1695         if (npdb.lt.10000) then                              
1696          write(liczba5,'(i1,i4)') 0,npdb
1697         else   
1698          write(liczba5,'(i5)') npdb
1699         endif
1700         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
1701       endif
1702       call pdbout(ee,titelloc1,ipdb)
1703       close(ipdb)
1704       return
1705       end
1706
1707 c-----------------------------------------------------------
1708       subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
1709       implicit real*8 (a-h,o-z)
1710       include 'DIMENSIONS'
1711 #ifdef MPI
1712       include 'mpif.h'
1713 #endif
1714       include 'COMMON.SBRIDGE'
1715       include 'COMMON.FFIELD'
1716       include 'COMMON.IOUNITS'
1717       include 'COMMON.DISTFIT'
1718       include 'COMMON.VAR'
1719       include 'COMMON.CHAIN'
1720       include 'COMMON.MINIM'
1721
1722       character*50 linia
1723       integer nf,ij(4)
1724       double precision var(maxvar),var2(maxvar)
1725       double precision time0,time1
1726       integer iff(maxres),ieval      
1727       double precision theta1(maxres),phi1(maxres),alph1(maxres),     
1728      &                 omeg1(maxres)                             
1729       
1730
1731        call var_to_geom(nvar,var)
1732        call chainbuild                                                           
1733        nhpb0=nhpb
1734        ind=0                                                                     
1735        do i=1,nres-3                                                             
1736          do j=i+3,nres                                                           
1737           ind=ind+1                                                              
1738           if ( iff(i).eq.1.and.iff(j).eq.1 ) then                                           
1739             d0(ind)=DIST(i,j)                                                     
1740             w(ind)=10.0                                                           
1741             nhpb=nhpb+1                                                           
1742             ihpb(nhpb)=i                                                          
1743             jhpb(nhpb)=j                                                          
1744             forcon(nhpb)=10.0                                                     
1745             dhpb(nhpb)=d0(ind)                                                    
1746           else
1747             w(ind)=0.0
1748           endif                                                                  
1749          enddo                                                                   
1750        enddo                                    
1751        call hpb_partition
1752
1753        do i=1,nres                                                               
1754         theta1(i)=theta(i)                                                      
1755         phi1(i)=phi(i)                                                          
1756         alph1(i)=alph(i)                                                        
1757         omeg1(i)=omeg(i)                                                        
1758        enddo                      
1759
1760        call var_to_geom(nvar,var2)
1761
1762        do i=1,nres                                                             
1763           if ( iff(i).eq.1 ) then                                           
1764               theta(i)=theta1(i)                                                      
1765               phi(i)=phi1(i)                                                          
1766               alph(i)=alph1(i)                                                        
1767               omeg(i)=omeg1(i)                       
1768           endif
1769        enddo
1770
1771        call chainbuild
1772 cd       call write_pdb(3,'combined structure',0d0)
1773 cd       time0=MPI_WTIME()
1774        
1775        NX=NRES-3                                                                 
1776        NY=((NRES-4)*(NRES-5))/2 
1777        call distfit(.true.,200)
1778    
1779 cd       time1=MPI_WTIME()
1780 cd       write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
1781
1782        ipot0=ipot
1783        maxmin0=maxmin
1784        maxfun0=maxfun
1785        wstrain0=wstrain
1786
1787        ipot=6
1788        maxmin=2000
1789        maxfun=5000
1790        call geom_to_var(nvar,var)
1791 cd       time0=MPI_WTIME()
1792        call minimize(etot,var,iretcode,nfun)                               
1793        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun   
1794
1795 cd       time1=MPI_WTIME()
1796 cd       write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
1797 cd     &         nfun/(time1-time0),' SOFT eval/s'
1798         call var_to_geom(nvar,var)
1799         call chainbuild
1800
1801
1802         iwsk=0
1803         nf=0
1804         if (iff(1).eq.1) then
1805           iwsk=1
1806           nf=nf+1
1807           ij(nf)=0
1808         endif
1809         do i=2,nres
1810            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
1811              iwsk=1
1812              nf=nf+1
1813              ij(nf)=i
1814            endif
1815            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
1816              iwsk=0
1817              nf=nf+1
1818              ij(nf)=i-1
1819            endif
1820         enddo
1821         if (iff(nres).eq.1) then
1822           nf=nf+1
1823           ij(nf)=nres
1824         endif
1825
1826
1827 cd        write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
1828 cd     &                     "select",ij(1),"-",ij(2),
1829 cd     &                         ",",ij(3),"-",ij(4)
1830 cd        call write_pdb(in_pdb,linia,etot)
1831
1832
1833        ipot=ipot0
1834        maxmin=maxmin0
1835        maxfun=maxfun0
1836 cd       time0=MPI_WTIME()
1837        call minimize(etot,var,iretcode,nfun)
1838 cd       write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
1839        ieval=nfun
1840
1841 cd       time1=MPI_WTIME()
1842 cd       write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
1843 cd     &         nfun/(time1-time0),' eval/s'
1844 cd       call var_to_geom(nvar,var)
1845 cd       call chainbuild
1846 cd       call write_pdb(6,'dist structure',etot)
1847
1848
1849        nhpb= nhpb0                                                                  
1850        link_start=1                                                            
1851        link_end=nhpb     
1852        wstrain=wstrain0
1853
1854        return
1855        end
1856 c-----------------------------------------------------------
1857       subroutine contact_cp(var,var2,iff,ieval,in_pdb)
1858       implicit real*8 (a-h,o-z)
1859       include 'DIMENSIONS'
1860       include 'COMMON.SBRIDGE'
1861       include 'COMMON.FFIELD'
1862       include 'COMMON.IOUNITS'
1863       include 'COMMON.DISTFIT'
1864       include 'COMMON.VAR'
1865       include 'COMMON.CHAIN'
1866       include 'COMMON.MINIM'
1867
1868       character*50 linia
1869       integer nf,ij(4)
1870       double precision energy(0:n_ene)
1871       double precision var(maxvar),var2(maxvar)
1872       double precision time0,time1
1873       integer iff(maxres),ieval      
1874       double precision theta1(maxres),phi1(maxres),alph1(maxres),     
1875      &                 omeg1(maxres)                             
1876       logical debug
1877       
1878       debug=.false.
1879 c      debug=.true.
1880       if (ieval.eq.-1) debug=.true.
1881
1882
1883 c
1884 c store selected dist. constrains from 1st structure
1885 c
1886 #ifdef OSF
1887 c     Intercept NaNs in the coordinates
1888 c      write(iout,*) (var(i),i=1,nvar)
1889       x_sum=0.D0
1890       do i=1,nvar
1891         x_sum=x_sum+var(i)
1892       enddo
1893       if (x_sum.ne.x_sum) then
1894         write(iout,*)" *** contact_cp : Found NaN in coordinates"
1895         call flush(iout) 
1896         print *," *** contact_cp : Found NaN in coordinates"
1897         return
1898       endif
1899 #endif
1900  
1901
1902        call var_to_geom(nvar,var)
1903        call chainbuild                                                           
1904        nhpb0=nhpb
1905        ind=0                                                                     
1906        do i=1,nres-3                                                             
1907          do j=i+3,nres                                                           
1908           ind=ind+1                                                              
1909           if ( iff(i).eq.1.and.iff(j).eq.1 ) then                                           
1910             d0(ind)=DIST(i,j)                                                     
1911             w(ind)=10.0                                                           
1912             nhpb=nhpb+1                                                           
1913             ihpb(nhpb)=i                                                          
1914             jhpb(nhpb)=j                                                          
1915             forcon(nhpb)=10.0                                                     
1916             dhpb(nhpb)=d0(ind)                                                    
1917           else
1918             w(ind)=0.0
1919           endif                                                                  
1920          enddo                                                                   
1921        enddo                                    
1922        call hpb_partition
1923
1924        do i=1,nres                                                               
1925         theta1(i)=theta(i)                                                      
1926         phi1(i)=phi(i)                                                          
1927         alph1(i)=alph(i)                                                        
1928         omeg1(i)=omeg(i)                                                        
1929        enddo                      
1930
1931 c
1932 c  freeze sec.elements from 2nd structure 
1933 c
1934        do i=1,nres
1935          mask_phi(i)=1
1936          mask_theta(i)=1
1937          mask_side(i)=1
1938        enddo
1939
1940        call var_to_geom(nvar,var2)
1941        call secondary2(debug)
1942        do j=1,nbfrag
1943         do i=bfrag(1,j),bfrag(2,j)
1944          mask(i)=0
1945          mask_phi(i)=0
1946          mask_theta(i)=0
1947         enddo
1948         if (bfrag(3,j).le.bfrag(4,j)) then 
1949          do i=bfrag(3,j),bfrag(4,j)
1950           mask(i)=0
1951           mask_phi(i)=0
1952           mask_theta(i)=0
1953          enddo
1954         else
1955          do i=bfrag(4,j),bfrag(3,j)
1956           mask(i)=0
1957           mask_phi(i)=0
1958           mask_theta(i)=0
1959          enddo
1960         endif
1961        enddo
1962        do j=1,nhfrag
1963         do i=hfrag(1,j),hfrag(2,j)
1964          mask(i)=0
1965          mask_phi(i)=0
1966          mask_theta(i)=0
1967         enddo
1968        enddo
1969        mask_r=.true.
1970
1971 c
1972 c      copy selected res from 1st to 2nd structure
1973 c
1974
1975        do i=1,nres                                                             
1976           if ( iff(i).eq.1 ) then                                           
1977               theta(i)=theta1(i)                                                      
1978               phi(i)=phi1(i)                                                          
1979               alph(i)=alph1(i)                                                        
1980               omeg(i)=omeg1(i)                       
1981           endif
1982        enddo
1983
1984       if(debug) then   
1985 c
1986 c     prepare description in linia variable
1987 c
1988         iwsk=0
1989         nf=0
1990         if (iff(1).eq.1) then
1991           iwsk=1
1992           nf=nf+1
1993           ij(nf)=1
1994         endif
1995         do i=2,nres
1996            if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
1997              iwsk=1
1998              nf=nf+1
1999              ij(nf)=i
2000            endif
2001            if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
2002              iwsk=0
2003              nf=nf+1
2004              ij(nf)=i-1
2005            endif
2006         enddo
2007         if (iff(nres).eq.1) then
2008           nf=nf+1
2009           ij(nf)=nres
2010         endif
2011
2012         write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
2013      &                     "SELECT",ij(1)-1,"-",ij(2)-1,
2014      &                         ",",ij(3)-1,"-",ij(4)-1
2015
2016       endif
2017 c
2018 c     run optimization
2019 c
2020       call contact_cp_min(var,ieval,in_pdb,linia,debug)
2021
2022       return
2023       end
2024
2025       subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
2026 c
2027 c    input : theta,phi,alph,omeg,in_pdb,linia,debug
2028 c    output : var,ieval
2029 c
2030       implicit real*8 (a-h,o-z)
2031       include 'DIMENSIONS'
2032 #ifdef MPI
2033       include 'mpif.h'
2034 #endif
2035       include 'COMMON.SBRIDGE'
2036       include 'COMMON.FFIELD'
2037       include 'COMMON.IOUNITS'
2038       include 'COMMON.DISTFIT'
2039       include 'COMMON.VAR'
2040       include 'COMMON.CHAIN'
2041       include 'COMMON.MINIM'
2042
2043       character*50 linia
2044       integer nf,ij(4)
2045       double precision energy(0:n_ene)
2046       double precision var(maxvar)
2047       double precision time0,time1
2048       integer ieval,info(3)      
2049       logical debug,fail,check_var,reduce,change
2050
2051        write(iout,'(a20,i6,a20)')
2052      &             '------------------',in_pdb,'-------------------'
2053
2054        if (debug) then
2055         call chainbuild
2056         call write_pdb(1000+in_pdb,'combined structure',0d0)
2057 #ifdef MPI
2058         time0=MPI_WTIME()
2059 #else
2060         time0=tcpu()
2061 #endif
2062        endif
2063        
2064 c
2065 c     run optimization of distances
2066 c     
2067 c     uses d0(),w() and mask() for frozen 2D
2068 c
2069 ctest---------------------------------------------
2070 ctest       NX=NRES-3                                                                 
2071 ctest       NY=((NRES-4)*(NRES-5))/2 
2072 ctest       call distfit(debug,5000)
2073
2074        do i=1,nres
2075          mask_side(i)=0
2076        enddo
2077  
2078        ipot01=ipot
2079        maxmin01=maxmin
2080        maxfun01=maxfun
2081 c       wstrain01=wstrain
2082        wsc01=wsc
2083        wscp01=wscp
2084        welec01=welec
2085        wvdwpp01=wvdwpp
2086 c      wang01=wang
2087        wscloc01=wscloc
2088        wtor01=wtor
2089        wtor_d01=wtor_d
2090
2091        ipot=6
2092        maxmin=2000
2093        maxfun=4000
2094 c       wstrain=1.0
2095        wsc=0.0
2096        wscp=0.0
2097        welec=0.0
2098        wvdwpp=0.0
2099 c      wang=0.0
2100        wscloc=0.0
2101        wtor=0.0
2102        wtor_d=0.0
2103
2104        call geom_to_var(nvar,var)
2105 cde       change=reduce(var)
2106        if (check_var(var,info)) then
2107           write(iout,*) 'cp_min error in input'
2108           print *,'cp_min error in input'
2109           return
2110        endif
2111
2112 cd       call etotal(energy(0))
2113 cd       call enerprint(energy(0))
2114 cd       call check_eint
2115
2116 #ifdef MPI
2117        time0=MPI_WTIME()
2118 #else
2119        time0=tcpu()
2120 #endif
2121 cdtest       call minimize(etot,var,iretcode,nfun)                               
2122 cdtest       write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun   
2123 #ifdef MPI
2124        time1=MPI_WTIME()
2125 #else
2126        time1=tcpu()
2127 #endif
2128
2129 cd       call etotal(energy(0))
2130 cd       call enerprint(energy(0))
2131 cd       call check_eint 
2132
2133        do i=1,nres
2134          mask_side(i)=1
2135        enddo
2136  
2137        ipot=ipot01
2138        maxmin=maxmin01
2139        maxfun=maxfun01
2140 c       wstrain=wstrain01
2141        wsc=wsc01
2142        wscp=wscp01
2143        welec=welec01
2144        wvdwpp=wvdwpp01
2145 c      wang=wang01
2146        wscloc=wscloc01
2147        wtor=wtor01
2148        wtor_d=wtor_d01
2149 ctest--------------------------------------------------
2150         
2151        if(debug) then
2152 #ifdef MPI
2153         time1=MPI_WTIME()
2154 #else
2155         time1=tcpu()
2156 #endif
2157         write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
2158         call write_pdb(2000+in_pdb,'distfit structure',0d0)
2159        endif
2160
2161
2162        ipot0=ipot
2163        maxmin0=maxmin
2164        maxfun0=maxfun
2165        wstrain0=wstrain
2166 c
2167 c      run soft pot. optimization 
2168 c         with constrains:
2169 c             nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition 
2170 c         and frozen 2D:
2171 c             mask_phi(),mask_theta(),mask_side(),mask_r
2172 c
2173        ipot=6
2174        maxmin=2000
2175        maxfun=4000
2176
2177 cde       change=reduce(var)
2178 cde       if (check_var(var,info)) write(iout,*) 'error before soft'
2179 #ifdef MPI
2180        time0=MPI_WTIME()
2181 #else
2182        time0=tcpu()
2183 #endif
2184        call minimize(etot,var,iretcode,nfun)                               
2185
2186        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun   
2187 #ifdef MPI
2188        time1=MPI_WTIME()
2189 #else
2190        time1=tcpu()
2191 #endif
2192        write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
2193      &         nfun/(time1-time0),' SOFT eval/s'
2194        if (debug) then
2195          call var_to_geom(nvar,var)
2196          call chainbuild
2197          call write_pdb(3000+in_pdb,'soft structure',etot)
2198        endif
2199 c
2200 c      run full UNRES optimization with constrains and frozen 2D
2201 c      the same variables as soft pot. optimizatio
2202 c
2203        ipot=ipot0
2204        maxmin=maxmin0
2205        maxfun=maxfun0
2206 c
2207 c check overlaps before calling full UNRES minim
2208 c
2209        call var_to_geom(nvar,var)
2210        call chainbuild
2211        call etotal(energy(0))
2212 #ifdef OSF
2213        write(iout,*) 'N7 ',energy(0)
2214        if (energy(0).ne.energy(0)) then
2215         write(iout,*) 'N7 error - gives NaN',energy(0)
2216        endif
2217 #endif
2218        ieval=1
2219        if (energy(1).eq.1.0d20) then
2220          write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
2221          call overlap_sc(fail)
2222          if(.not.fail) then
2223            call etotal(energy(0))
2224            ieval=ieval+1
2225            write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
2226          else
2227            mask_r=.false.
2228            nhpb= nhpb0
2229            link_start=1
2230            link_end=nhpb
2231            wstrain=wstrain0
2232            return
2233          endif
2234        endif
2235        call flush(iout)
2236 c
2237 cdte       time0=MPI_WTIME()
2238 cde       change=reduce(var)
2239 cde       if (check_var(var,info)) then 
2240 cde         write(iout,*) 'error before mask dist'
2241 cde         call var_to_geom(nvar,var)
2242 cde         call chainbuild
2243 cde         call write_pdb(10000+in_pdb,'before mask dist',etot)
2244 cde       endif
2245 cdte       call minimize(etot,var,iretcode,nfun)
2246 cdte       write(iout,*)'SUMSL MASK DIST return code is',iretcode,
2247 cdte     &                          ' eval ',nfun
2248 cdte       ieval=ieval+nfun
2249 cdte
2250 cdte       time1=MPI_WTIME()
2251 cdte       write (iout,'(a,f6.2,f8.2,a)') 
2252 cdte     &        ' Time for mask dist min.',time1-time0,
2253 cdte     &         nfun/(time1-time0),'  eval/s'
2254 cdte       call flush(iout)
2255        if (debug) then
2256          call var_to_geom(nvar,var)
2257          call chainbuild
2258          call write_pdb(4000+in_pdb,'mask dist',etot)
2259        endif
2260 c
2261 c      switch off freezing of 2D and 
2262 c      run full UNRES optimization with constrains 
2263 c
2264        mask_r=.false.
2265 #ifdef MPI
2266        time0=MPI_WTIME()
2267 #else
2268        time0=tcpu()
2269 #endif
2270 cde       change=reduce(var)
2271 cde       if (check_var(var,info)) then 
2272 cde         write(iout,*) 'error before dist'
2273 cde         call var_to_geom(nvar,var)
2274 cde         call chainbuild
2275 cde         call write_pdb(11000+in_pdb,'before dist',etot)
2276 cde       endif
2277
2278        call minimize(etot,var,iretcode,nfun)
2279
2280 cde        change=reduce(var)
2281 cde        if (check_var(var,info)) then 
2282 cde          write(iout,*) 'error after dist',ico
2283 cde          call var_to_geom(nvar,var)
2284 cde          call chainbuild
2285 cde          call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
2286 cde        endif
2287        write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
2288        ieval=ieval+nfun
2289
2290 #ifdef MPI
2291        time1=MPI_WTIME()
2292 #else
2293        time1=tcpu()
2294 #endif
2295        write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,
2296      &         nfun/(time1-time0),'  eval/s'
2297 cde       call etotal(energy(0))
2298 cde       write(iout,*) 'N7 after dist',energy(0)
2299        call flush(iout)
2300
2301        if (debug) then
2302         call var_to_geom(nvar,var)
2303         call chainbuild
2304         call write_pdb(in_pdb,linia,etot)
2305        endif
2306 c
2307 c      reset constrains
2308 c
2309        nhpb= nhpb0                                                                 
2310        link_start=1                                                            
2311        link_end=nhpb     
2312        wstrain=wstrain0
2313
2314        return
2315        end
2316 c--------------------------------------------------------
2317       subroutine softreg
2318       implicit real*8 (a-h,o-z)
2319       include 'DIMENSIONS'
2320 #ifdef MPI
2321       include 'mpif.h'
2322 #endif
2323       include 'COMMON.GEO'
2324       include 'COMMON.CHAIN'
2325       include 'COMMON.IOUNITS'
2326       include 'COMMON.VAR'
2327       include 'COMMON.CONTROL'
2328       include 'COMMON.SBRIDGE'
2329       include 'COMMON.FFIELD'
2330       include 'COMMON.MINIM'
2331       include 'COMMON.INTERACT'
2332 c
2333       include 'COMMON.DISTFIT'       
2334       integer iff(maxres)
2335       double precision time0,time1
2336       double precision energy(0:n_ene),ee
2337       double precision var(maxvar)
2338       integer ieval
2339 c
2340       logical debug,ltest,fail
2341       character*50 linia
2342 c
2343       linia='test'
2344       debug=.true.
2345       in_pdb=0
2346
2347
2348
2349 c------------------------
2350 c
2351 c  freeze sec.elements 
2352 c
2353        do i=1,nres
2354          mask_phi(i)=1
2355          mask_theta(i)=1
2356          mask_side(i)=1
2357          iff(i)=0
2358        enddo
2359
2360        do j=1,nbfrag
2361         do i=bfrag(1,j),bfrag(2,j)
2362          mask_phi(i)=0
2363          mask_theta(i)=0
2364          iff(i)=1
2365         enddo
2366         if (bfrag(3,j).le.bfrag(4,j)) then 
2367          do i=bfrag(3,j),bfrag(4,j)
2368           mask_phi(i)=0
2369           mask_theta(i)=0
2370           iff(i)=1
2371          enddo
2372         else
2373          do i=bfrag(4,j),bfrag(3,j)
2374           mask_phi(i)=0
2375           mask_theta(i)=0
2376           iff(i)=1
2377          enddo
2378         endif
2379        enddo
2380        do j=1,nhfrag
2381         do i=hfrag(1,j),hfrag(2,j)
2382          mask_phi(i)=0
2383          mask_theta(i)=0
2384          iff(i)=1
2385         enddo
2386        enddo
2387        mask_r=.true.
2388
2389
2390
2391        nhpb0=nhpb
2392 c
2393 c store dist. constrains
2394 c
2395        do i=1,nres-3                                                             
2396          do j=i+3,nres                                                           
2397            if ( iff(i).eq.1.and.iff(j).eq.1 ) then
2398             nhpb=nhpb+1                                                           
2399             ihpb(nhpb)=i                                                          
2400             jhpb(nhpb)=j                                                          
2401             forcon(nhpb)=0.1                                                     
2402             dhpb(nhpb)=DIST(i,j)
2403            endif
2404          enddo                                                                   
2405        enddo                                    
2406        call hpb_partition
2407
2408        if (debug) then
2409         call chainbuild
2410         call write_pdb(100+in_pdb,'input reg. structure',0d0)
2411        endif
2412        
2413
2414        ipot0=ipot
2415        maxmin0=maxmin
2416        maxfun0=maxfun
2417        wstrain0=wstrain
2418        wang0=wang
2419 c
2420 c      run soft pot. optimization 
2421 c
2422        ipot=6
2423        wang=3.0
2424        maxmin=2000
2425        maxfun=4000
2426        call geom_to_var(nvar,var)
2427 #ifdef MPI
2428        time0=MPI_WTIME()
2429 #else
2430        time0=tcpu()
2431 #endif
2432        call minimize(etot,var,iretcode,nfun)                               
2433
2434        write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun   
2435 #ifdef MPI
2436        time1=MPI_WTIME()
2437 #else
2438        time1=tcpu()
2439 #endif
2440        write (iout,'(a,f6.2,f8.2,a)')'  Time for soft min.',time1-time0,
2441      &         nfun/(time1-time0),' SOFT eval/s'
2442        if (debug) then
2443          call var_to_geom(nvar,var)
2444          call chainbuild
2445          call write_pdb(300+in_pdb,'soft structure',etot)
2446        endif
2447 c
2448 c      run full UNRES optimization with constrains and frozen 2D
2449 c      the same variables as soft pot. optimizatio
2450 c
2451        ipot=ipot0
2452        wang=wang0
2453        maxmin=maxmin0
2454        maxfun=maxfun0
2455 #ifdef MPI
2456        time0=MPI_WTIME()
2457 #else
2458        time0=tcpu()
2459 #endif
2460        call minimize(etot,var,iretcode,nfun)
2461        write(iout,*)'SUMSL MASK DIST return code is',iretcode,
2462      &                          ' eval ',nfun
2463        ieval=nfun
2464
2465 #ifdef MPI
2466        time1=MPI_WTIME()
2467 #else
2468        time1=tcpu()
2469 #endif
2470        write (iout,'(a,f6.2,f8.2,a)') 
2471      &        '  Time for mask dist min.',time1-time0,
2472      &         nfun/(time1-time0),'  eval/s'
2473        if (debug) then
2474          call var_to_geom(nvar,var)
2475          call chainbuild
2476          call write_pdb(400+in_pdb,'mask & dist',etot)
2477        endif
2478 c
2479 c      switch off constrains and 
2480 c      run full UNRES optimization with frozen 2D 
2481 c
2482
2483 c
2484 c      reset constrains
2485 c
2486        nhpb_c=nhpb
2487        nhpb=nhpb0                                                                  
2488        link_start=1                                                            
2489        link_end=nhpb     
2490        wstrain=wstrain0
2491
2492
2493 #ifdef MPI
2494        time0=MPI_WTIME()
2495 #else
2496        time0=tcpu()
2497 #endif
2498        call minimize(etot,var,iretcode,nfun)
2499        write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
2500        ieval=ieval+nfun
2501
2502 #ifdef MPI
2503        time1=MPI_WTIME()
2504 #else
2505        time1=tcpu()
2506 #endif
2507        write (iout,'(a,f6.2,f8.2,a)')'  Time for mask min.',time1-time0,
2508      &         nfun/(time1-time0),'  eval/s'
2509
2510
2511        if (debug) then
2512         call var_to_geom(nvar,var)
2513         call chainbuild
2514         call write_pdb(500+in_pdb,'mask 2d frozen',etot)
2515        endif
2516
2517        mask_r=.false.
2518
2519
2520 c
2521 c      run full UNRES optimization with constrains and NO frozen 2D
2522 c
2523
2524        nhpb=nhpb_c                                                                  
2525        link_start=1                                                            
2526        link_end=nhpb     
2527        maxfun=maxfun0/5
2528
2529        do ico=1,5
2530
2531        wstrain=wstrain0/ico
2532
2533 #ifdef MPI
2534        time0=MPI_WTIME()
2535 #else
2536        time0=tcpu()
2537 #endif
2538        call minimize(etot,var,iretcode,nfun)
2539        write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2540      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2541      &                          ' eval ',nfun
2542        ieval=nfun
2543
2544 #ifdef MPI
2545        time1=MPI_WTIME()
2546 #else
2547        time0=tcpu()
2548 #endif
2549        write (iout,'(a,f6.2,f8.2,a)') 
2550      &        '  Time for dist min.',time1-time0,
2551      &         nfun/(time1-time0),'  eval/s'
2552        if (debug) then
2553          call var_to_geom(nvar,var)
2554          call chainbuild
2555          call write_pdb(600+in_pdb+ico,'dist cons',etot)
2556        endif
2557
2558        enddo
2559 c
2560        nhpb=nhpb0                                                                  
2561        link_start=1                                                            
2562        link_end=nhpb     
2563        wstrain=wstrain0
2564        maxfun=maxfun0
2565
2566
2567 c
2568       if (minim) then
2569 #ifdef MPI
2570        time0=MPI_WTIME()
2571 #else
2572        time0=tcpu()
2573 #endif
2574        call minimize(etot,var,iretcode,nfun)
2575        write(iout,*)'------------------------------------------------'
2576        write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
2577      &  '+ DIST eval',ieval
2578       
2579 #ifdef MPI
2580        time1=MPI_WTIME()
2581 #else
2582        time1=tcpu()
2583 #endif
2584        write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
2585      &         nfun/(time1-time0),' eval/s'
2586
2587
2588        call var_to_geom(nvar,var)
2589        call chainbuild        
2590        call write_pdb(999,'full min',etot)
2591       endif
2592        
2593       return
2594       end
2595
2596
2597       subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
2598       implicit real*8 (a-h,o-z)
2599       include 'DIMENSIONS'
2600 #ifdef MPI
2601       include 'mpif.h'
2602 #endif
2603       include 'COMMON.GEO'
2604       include 'COMMON.VAR'
2605       include 'COMMON.INTERACT'
2606       include 'COMMON.IOUNITS'
2607       include 'COMMON.DISTFIT'
2608       include 'COMMON.SBRIDGE'
2609       include 'COMMON.CONTROL'
2610       include 'COMMON.FFIELD'
2611       include 'COMMON.MINIM'
2612       include 'COMMON.CHAIN'
2613       double precision time0,time1
2614       double precision energy(0:n_ene),ee
2615       double precision var(maxvar)
2616       integer jdata(5),isec(maxres)
2617 c
2618       jdata(1)=i1
2619       jdata(2)=i2
2620       jdata(3)=i3
2621       jdata(4)=i4
2622       jdata(5)=i5
2623
2624       call secondary2(.false.)
2625
2626       do i=1,nres
2627           isec(i)=0
2628       enddo
2629       do j=1,nbfrag
2630        do i=bfrag(1,j),bfrag(2,j)
2631           isec(i)=1
2632        enddo
2633        do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
2634           isec(i)=1
2635        enddo
2636       enddo
2637       do j=1,nhfrag
2638        do i=hfrag(1,j),hfrag(2,j)
2639           isec(i)=2
2640        enddo
2641       enddo
2642
2643 c
2644 c cut strands at the ends
2645 c
2646       if (jdata(2)-jdata(1).gt.3) then
2647        jdata(1)=jdata(1)+1
2648        jdata(2)=jdata(2)-1
2649        if (jdata(3).lt.jdata(4)) then
2650           jdata(3)=jdata(3)+1
2651           jdata(4)=jdata(4)-1
2652        else
2653           jdata(3)=jdata(3)-1
2654           jdata(4)=jdata(4)+1    
2655        endif
2656       endif
2657
2658 cv      call chainbuild
2659 cv      call etotal(energy(0))
2660 cv      etot=energy(0)
2661 cv      write(iout,*) nnt,nct,etot
2662 cv      call write_pdb(ij*100,'first structure',etot)
2663 cv      write(iout,*) 'N16 test',(jdata(i),i=1,5)
2664
2665 c------------------------
2666 c      generate constrains 
2667 c
2668        ishift=jdata(5)-2
2669        if(ishift.eq.0) ishift=-2
2670        nhpb0=nhpb
2671        call chainbuild                                                           
2672        do i=jdata(1),jdata(2)                                                             
2673         isec(i)=-1
2674         if(jdata(4).gt.jdata(3))then
2675          do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
2676             isec(j)=-1
2677 cd            print *,i,j,j+ishift
2678             nhpb=nhpb+1                                                           
2679             ihpb(nhpb)=i                                                          
2680             jhpb(nhpb)=j                                                          
2681             forcon(nhpb)=1000.0                                                     
2682             dhpb(nhpb)=DIST(i,j+ishift)
2683          enddo               
2684         else
2685          do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
2686             isec(j)=-1
2687 cd            print *,i,j,j+ishift
2688             nhpb=nhpb+1                                                           
2689             ihpb(nhpb)=i                                                          
2690             jhpb(nhpb)=j                                                          
2691             forcon(nhpb)=1000.0                                                     
2692             dhpb(nhpb)=DIST(i,j+ishift)
2693          enddo
2694         endif                                                    
2695        enddo      
2696
2697        do i=nnt,nct-2
2698          do j=i+2,nct
2699            if(isec(i).gt.0.or.isec(j).gt.0) then
2700 cd            print *,i,j
2701             nhpb=nhpb+1
2702             ihpb(nhpb)=i
2703             jhpb(nhpb)=j
2704             forcon(nhpb)=0.1
2705             dhpb(nhpb)=DIST(i,j)
2706            endif
2707          enddo
2708        enddo
2709                               
2710        call hpb_partition
2711
2712        call geom_to_var(nvar,var)       
2713        maxfun0=maxfun
2714        wstrain0=wstrain
2715        maxfun=4000/5
2716
2717        do ico=1,5
2718
2719         wstrain=wstrain0/ico
2720
2721 cv        time0=MPI_WTIME()
2722         call minimize(etot,var,iretcode,nfun)
2723         write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2724      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2725      &                          ' eval ',nfun
2726         ieval=ieval+nfun
2727 cv        time1=MPI_WTIME()
2728 cv       write (iout,'(a,f6.2,f8.2,a)') 
2729 cv     &        '  Time for dist min.',time1-time0,
2730 cv     &         nfun/(time1-time0),'  eval/s'
2731 cv         call var_to_geom(nvar,var)
2732 cv         call chainbuild
2733 cv         call write_pdb(ij*100+ico,'dist cons',etot)
2734
2735        enddo
2736 c
2737        nhpb=nhpb0                                                                  
2738        call hpb_partition
2739        wstrain=wstrain0
2740        maxfun=maxfun0
2741 c
2742 cd      print *,etot
2743       wscloc0=wscloc
2744       wscloc=10.0
2745       call sc_move(nnt,nct,100,100d0,nft_sc,etot)
2746       wscloc=wscloc0
2747 cv      call chainbuild
2748 cv      call etotal(energy(0))
2749 cv      etot=energy(0)
2750 cv      call write_pdb(ij*100+10,'sc_move',etot)
2751 cd      call intout
2752 cd      print *,nft_sc,etot
2753
2754       return
2755       end
2756
2757       subroutine beta_zip(i1,i2,ieval,ij)
2758       implicit real*8 (a-h,o-z)
2759       include 'DIMENSIONS'
2760 #ifdef MPI
2761       include 'mpif.h'
2762 #endif
2763       include 'COMMON.GEO'
2764       include 'COMMON.VAR'
2765       include 'COMMON.INTERACT'
2766       include 'COMMON.IOUNITS'
2767       include 'COMMON.DISTFIT'
2768       include 'COMMON.SBRIDGE'
2769       include 'COMMON.CONTROL'
2770       include 'COMMON.FFIELD'
2771       include 'COMMON.MINIM'
2772       include 'COMMON.CHAIN'
2773       double precision time0,time1
2774       double precision energy(0:n_ene),ee
2775       double precision var(maxvar)
2776       character*10 test
2777
2778 cv      call chainbuild
2779 cv      call etotal(energy(0))
2780 cv      etot=energy(0)
2781 cv      write(test,'(2i5)') i1,i2
2782 cv      call write_pdb(ij*100,test,etot)
2783 cv      write(iout,*) 'N17 test',i1,i2,etot,ij
2784
2785 c
2786 c      generate constrains 
2787 c
2788        nhpb0=nhpb
2789        nhpb=nhpb+1                                                           
2790        ihpb(nhpb)=i1                                                          
2791        jhpb(nhpb)=i2                                                          
2792        forcon(nhpb)=1000.0                                                     
2793        dhpb(nhpb)=4.0
2794                               
2795        call hpb_partition
2796
2797        call geom_to_var(nvar,var)       
2798        maxfun0=maxfun
2799        wstrain0=wstrain
2800        maxfun=1000/5
2801
2802        do ico=1,5
2803         wstrain=wstrain0/ico
2804 cv        time0=MPI_WTIME()
2805         call minimize(etot,var,iretcode,nfun)
2806         write(iout,'(a10,f6.3,a14,i3,a6,i5)')
2807      &   ' SUMSL DIST',wstrain,' return code is',iretcode,
2808      &                          ' eval ',nfun
2809         ieval=ieval+nfun
2810 cv        time1=MPI_WTIME()
2811 cv       write (iout,'(a,f6.2,f8.2,a)') 
2812 cv     &        '  Time for dist min.',time1-time0,
2813 cv     &         nfun/(time1-time0),'  eval/s'
2814 c do not comment the next line
2815          call var_to_geom(nvar,var)
2816 cv         call chainbuild
2817 cv         call write_pdb(ij*100+ico,'dist cons',etot)
2818        enddo
2819
2820        nhpb=nhpb0                                                                  
2821        call hpb_partition
2822        wstrain=wstrain0
2823        maxfun=maxfun0
2824
2825 cv      call etotal(energy(0))
2826 cv      etot=energy(0)
2827 cv      write(iout,*) 'N17 test end',i1,i2,etot,ij
2828
2829
2830       return
2831       end