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