multichain wham ignores conf. energies>10E6 and speed up of starts
[unres.git] / source / wham / src / enecalc1.F
1       subroutine enecalc(islice,*)
2       implicit none
3       include "DIMENSIONS"
4       include "DIMENSIONS.ZSCOPT"
5       include "DIMENSIONS.FREE"
6 #ifdef MPI
7       include "mpif.h"
8       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
9       include "COMMON.MPI"
10 #endif
11       include "COMMON.CHAIN"
12       include "COMMON.IOUNITS"
13       include "COMMON.PROTFILES"
14       include "COMMON.NAMES"
15       include "COMMON.VAR"
16       include "COMMON.SBRIDGE"
17       include "COMMON.GEO"
18       include "COMMON.FFIELD"
19       include "COMMON.ENEPS"
20       include "COMMON.LOCAL"
21       include "COMMON.WEIGHTS"
22       include "COMMON.INTERACT"
23       include "COMMON.FREE"
24       include "COMMON.ENERGIES"
25       include "COMMON.CONTROL"
26       include "COMMON.TORCNSTR"
27       character*64 nazwa
28       character*80 bxname
29       character*3 liczba
30       double precision qwolynes
31       external qwolynes
32       integer errmsg_count,maxerrmsg_count /100/ 
33       double precision rmsnat,gyrate
34       external rmsnat,gyrate
35       double precision tole /1.0d-1/
36       integer i,itj,ii,iii,j,k,l,licz,scme,itmp
37       integer ires
38       integer ir,ib,ipar,iparm
39       integer iscor,islice
40       real*4 csingle(3,maxres2)
41       double precision energ
42       integer ilen,iroof
43       external ilen,iroof
44       double precision energia(0:max_ene),rmsdev,efree,eini
45       double precision fT(6),quot,quotl,kfacl,kfac /2.4d0/,T0 /3.0d2/
46       double precision tt
47       integer snk_p(MaxR,MaxT_h,Max_parm)
48       logical lerr
49       character*128 bprotfile_temp
50       integer scount_t(0:maxprocs-1)
51       call opentmp(islice,ientout,bprotfile_temp)
52       iii=0
53       ii=0
54       errmsg_count=0
55       write (iout,*) "enecalc: nparmset ",nparmset
56 #ifdef MPI
57       do iparm=1,nParmSet
58         do ib=1,nT_h(iparm)
59           do i=1,nR(ib,iparm)
60             snk_p(i,ib,iparm)=0
61           enddo
62         enddo
63       enddo
64       do i=indstart(me1),indend(me1)
65 #else
66       do iparm=1,nParmSet
67         do ib=1,nT_h(iparm)
68           do i=1,nR(ib,iparm)
69             snk(i,ib,iparm)=0
70           enddo
71         enddo
72       enddo
73       do i=1,ntot
74 #endif
75         read(ientout,rec=i,err=101) 
76      &    ((csingle(l,k),l=1,3),k=1,nres),
77      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
78      &    nss,(ihpb(k),jhpb(k),k=1,nss),
79      &    eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
80 #ifdef DEBUG
81         write (iout,*) "homol_nset",homol_nset,
82      &    "  i",i," iR",iR," ib",ib," iset",iset
83 #endif
84         if (homol_nset.gt.1) iset=iR
85 #ifdef DEBUG
86         write (iout,*) "homol_nset",homol_nset,
87      &    "  i",i," iR",iR," ib",ib," iset",iset
88 #endif
89 cc       write(iout,*), 'NAWEJ',i,eini
90          if (indpdb.gt.0) then
91            do k=1,nres
92              do l=1,3
93                c(l,k)=csingle(l,k)
94              enddo
95            enddo
96            do k=nnt,nct
97              do l=1,3
98                c(l,k+nres)=csingle(l,k+nres)
99              enddo
100            enddo
101            q(nQ+1,iii+1)=rmsnat(iii+1)
102          endif
103          q(nQ+2,iii+1)=gyrate(iii+1)
104 c        fT=T0*beta_h(ib,ipar)*1.987D-3
105 c        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
106         if (rescale_mode.eq.1) then
107           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
108 #if defined(FUNCTH)
109           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
110           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
111 #elif defined(FUNCT)
112           ft(6)=quot
113 #else
114           ft(6)=1.0d0
115 #endif
116           quotl=1.0d0
117           kfacl=1.0d0
118           do l=1,5
119             quotl=quotl*quot
120             kfacl=kfacl*kfac
121             fT(l)=kfacl/(kfacl-1.0d0+quotl)
122           enddo
123         else if (rescale_mode.eq.2) then
124           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
125 #if defined(FUNCTH)
126           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
127           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
128 #elif defined(FUNCT)
129           ft(6)=quot
130 #else
131           ft(6)=1.0d0
132 #endif
133           quotl=1.0d0
134           do l=1,5
135             quotl=quotl*quot
136             fT(l)=1.12692801104297249644d0/
137      &         dlog(dexp(quotl)+dexp(-quotl))
138           enddo
139         else if (rescale_mode.eq.0) then
140           do l=1,5
141             fT(l)=1.0d0
142           enddo
143         else
144           write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",
145      &     rescale_mode
146           call flush(iout)
147           return1
148         endif
149
150 c        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
151 c     &   " kfac",kfac,"quot",quot," fT",fT
152         do j=1,2*nres
153           do k=1,3
154             c(k,j)=csingle(k,j)
155           enddo
156         enddo
157         call int_from_cart1(.false.)
158         ii=ii+1
159         do iparm=1,nparmset
160
161         call restore_parm(iparm)
162 #ifdef DEBUG
163             write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
164      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
165      &      wtor_d,wsccor,wbond
166 #endif
167         call etotal(energia(0),fT)
168         if (constr_homology.gt.0) energia(0)=energia(0)+
169      &   waga_homology(iset)*energia(22)
170 c        write (iout,*) "constr_homology",constr_homology," iset",iset,
171 c     &   " waga_homology",waga_homology(iset)
172 #ifdef DEBUG
173         write (iout,*) "Conformation",i
174         call enerprint(energia(0),fT)
175 c        write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
176 c        write (iout,*) "ftors",ftors
177 c      write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
178 c     & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
179 c      do ires=1,nres
180 c        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
181 c     &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
182 c     &    (c(j,ires+nres),j=1,3)
183 c      enddo
184         call intout
185 #endif
186         if (energia(0).ge.1.0d20) then
187           write (iout,*) "NaNs detected in some of the energy",
188      &     " components for conformation",ii+1
189           write (iout,*) "The Cartesian geometry is:"
190           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
191           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
192           write (iout,*) "The internal geometry is:"
193 c          call intout
194 c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
195           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
196           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
197           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
198           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
199           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
200           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
201           write (iout,*) "The components of the energy are:"
202           call enerprint(energia(0),fT)
203           write (iout,*) 
204      &      "This conformation WILL NOT be added to the database."
205           call flush(iout)
206           goto 121
207         else 
208 #ifdef DEBUG
209           if (ipar.eq.iparm) write (iout,*) i,iparm,
210      &      1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
211 #endif
212           if (ipar.eq.iparm .and. einicheck.gt.0 .and. 
213      &      dabs(eini-energia(0)).gt.tole) then
214             if (errmsg_count.le.maxerrmsg_count) then
215               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)') 
216      &         "Warning: energy differs remarkably from ",
217      &         " the value read in: ",energia(0),eini," point",
218      &         iii+1,indstart(me1)+iii," T",
219      &         1.0d0/(1.987D-3*beta_h(ib,ipar))
220 #ifdef DEBUG
221              call enerprint(energia(0),fT)
222 #endif
223 #ifdef DEBUG
224              write (iout,'(4f10.5,2i5)') 0.0,energia(0),0.0,
225      &       1.0d0/(beta_h(ib,ipar)*1.987D-3),
226      &       0,0
227              write(iout,'(8f10.5)')
228      &       ((c(l,k),l=1,3),k=1,nres),
229      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
230              itmp=ipdb
231              ipdb=iout
232              call pdbout(iii+1,beta_h(ib,ipar),
233      &                   eini,energia(0),0.0d0,rmsdev)
234              write (iout,*)
235              ipdb=itmp
236 #endif
237               errmsg_count=errmsg_count+1
238               if (errmsg_count.gt.maxerrmsg_count) 
239      &          write (iout,*) "Too many warning messages"
240               if (einicheck.gt.1) then
241                 write (iout,*) "Calculation stopped."
242                 call flush(iout)
243 #ifdef MPI
244                 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
245 #endif
246                 call flush(iout)
247                 return1
248               endif
249             endif
250           endif
251           potE(iii+1,iparm)=energia(0)
252           do k=1,max_ene
253             enetb(k,iii+1,iparm)=energia(k)
254           enddo
255 c          write(iout,*) "iCHUJ TU STRZELI",i,iii,entfac(i)
256 c          call enerprint(energia(0),fT)
257 #ifdef DEBUG
258           write (iout,'(2i5,f10.1,3e15.5)') i,iii,
259      &     1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
260           call enerprint(energia(0),fT)
261           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
262           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
263           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
264           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
265           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
266           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
267           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
268           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
269           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
270           write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
271           write (iout,'(f10.5,i10)') rmsdev,iscor
272           call enerprint(energia(0),fT)
273         write(liczba,'(bz,i3.3)') me
274         nazwa="test"//liczba//".pdb"
275         write (iout,*) "pdb file",nazwa
276         open (ipdb,file=nazwa,position="append")
277         call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
278         close(ipdb)
279 #endif
280         endif
281
282         enddo ! iparm
283
284         iii=iii+1
285         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0) q(1,iii)=qwolynes(0,0)
286         write (ientout,rec=iii) 
287      &   ((csingle(l,k),l=1,3),k=1,nres),
288      &   ((csingle(l,k+nres),l=1,3),k=nnt,nct),
289      &   nss,(ihpb(k),jhpb(k),k=1,nss),
290      &   potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
291 c        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
292 #ifdef MPI
293         if (separate_parset) then
294           snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
295         else
296           snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
297         endif
298 c        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
299 c     &   " snk",snk_p(iR,ib,ipar)
300 #else
301         snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
302 #endif
303   121   continue
304       enddo   
305 #ifdef MPI
306       scount(me)=iii 
307       write (iout,*) "Me",me," scount",scount(me)
308       call flush(iout)
309 c  Master gathers updated numbers of conformations written by all procs.
310       scme = scount(me)
311       call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_t(0), 1, 
312      &  MPI_INTEGER, WHAM_COMM, IERROR)
313       do k=0,nprocs-1
314         scount(k) = scount_t(k)
315       enddo
316       indstart(0)=1
317       indend(0)=scount(0)
318       do i=1, Nprocs-1
319         indstart(i)=indend(i-1)+1
320         indend(i)=indstart(i)+scount(i)-1
321       enddo
322       write (iout,*)
323       write (iout,*) "Revised conformation counts"
324       do i=0,nprocs1-1
325         write (iout,'(a,i5,a,i7,a,i7,a,i7)')
326      &    "Processor",i," indstart",indstart(i),
327      &    " indend",indend(i)," count",scount(i)
328       enddo
329       call flush(iout)
330       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),
331      &  MaxR*MaxT_h*nParmSet,
332      &  MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
333 #endif
334       stot(islice)=0
335       do iparm=1,nParmSet
336         do ib=1,nT_h(iparm)
337           do i=1,nR(ib,iparm)
338             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
339           enddo
340         enddo
341       enddo
342       write (iout,*) "Revised SNK"
343       do iparm=1,nParmSet
344         do ib=1,nT_h(iparm)
345           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)') 
346      &     iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),
347      &     (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
348           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
349         enddo
350       enddo
351       write (iout,'("Total",i10)') stot(islice)
352       call flush(iout)
353       return
354   101 write (iout,*) "Error in scratchfile."
355       call flush(iout)
356       return1
357       end
358 c------------------------------------------------------------------------------
359       subroutine write_dbase(islice,*)
360       implicit none
361       include "DIMENSIONS"
362       include "DIMENSIONS.ZSCOPT"
363       include "DIMENSIONS.FREE"
364       include "DIMENSIONS.COMPAR"
365 #ifdef MPI
366       include "mpif.h"
367       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
368       include "COMMON.MPI"
369 #endif
370       include "COMMON.CONTROL"
371       include "COMMON.CHAIN"
372       include "COMMON.IOUNITS"
373       include "COMMON.PROTFILES"
374       include "COMMON.NAMES"
375       include "COMMON.VAR"
376       include "COMMON.SBRIDGE"
377       include "COMMON.GEO"
378       include "COMMON.FFIELD"
379       include "COMMON.ENEPS"
380       include "COMMON.LOCAL"
381       include "COMMON.WEIGHTS"
382       include "COMMON.INTERACT"
383       include "COMMON.FREE"
384       include "COMMON.ENERGIES"
385       include "COMMON.COMPAR"
386       include "COMMON.PROT"
387       character*64 nazwa
388       character*80 bxname,cxname
389       character*128 bprotfile_temp
390       character*3 liczba,licz
391       character*2 licz2
392       integer i,itj,ii,iii,j,k,l
393       integer ixdrf,iret
394       integer iscor,islice
395       double precision rmsdev,efree,eini
396       real*4 csingle(3,maxres2)
397       double precision energ
398       integer ilen,iroof
399       external ilen,iroof
400       integer ir,ib,iparm
401       write (licz2,'(bz,i2.2)') islice
402       call opentmp(islice,ientout,bprotfile_temp)
403       write (iout,*) "bprotfile_temp ",bprotfile_temp
404       call flush(iout)
405       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0 
406      &   .and. ensembles.eq.0) then
407         close(ientout,status="delete")
408         return
409       endif
410 #ifdef MPI
411       write (liczba,'(bz,i3.3)') me
412       if (bxfile .or. cxfile .or. ensembles.gt.0) then
413         if (.not.separate_parset) then
414           bxname = prefix(:ilen(prefix))//liczba//".bx"
415         else
416           write (licz,'(bz,i3.3)') myparm
417           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
418         endif
419         open (ientin,file=bxname,status="unknown",
420      &    form="unformatted",access="direct",recl=lenrec1)
421       endif
422 #else
423       if (bxfile .or. cxfile .or. ensembles.gt.0) then
424         if (nslice.eq.1) then
425           bxname = prefix(:ilen(prefix))//".bx"
426         else
427           bxname = prefix(:ilen(prefix))//
428      &           "_slice_"//licz2//".bx"
429         endif
430         open (ientin,file=bxname,status="unknown",
431      &    form="unformatted",access="direct",recl=lenrec1)
432         write (iout,*) "Calculating energies; writing geometry",
433      & " and energy components to ",bxname(:ilen(bxname))
434       endif
435 #if (defined(AIX) && !defined(JUBL))
436         call xdrfopen_(ixdrf,cxname, "w", iret)
437 #else
438         call xdrfopen(ixdrf,cxname, "w", iret)
439 #endif
440         if (iret.eq.0) then
441           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
442           cxfile=.fale.
443         endif
444       endif 
445 #endif
446       if (indpdb.gt.0) then
447         if (nslice.eq.1) then
448 #ifdef MPI
449          if (.not.separate_parset) then
450            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
451      &       //liczba//'.stat'
452          else
453            write (licz,'(bz,i3.3)') myparm
454            statname=prefix(:ilen(prefix))//'_par'//licz//'_'//
455      &      pot(:ilen(pot))//liczba//'.stat'
456          endif
457
458 #else
459           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
460 #endif
461         else
462 #ifdef MPI
463          if (.not.separate_parset) then
464           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
465      &      "_slice_"//licz2//liczba//'.stat'
466          else
467           write (licz,'(bz,i3.3)') myparm
468           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
469      &      '_par'//licz//"_slice_"//licz2//liczba//'.stat'
470          endif
471 #else
472           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
473      &      //"_slice_"//licz2//'.stat'
474 #endif
475         endif
476         open(istat,file=statname,status="unknown")
477       endif
478
479 #ifdef MPI
480       do i=1,scount(me)
481 #else
482       do i=1,ntot(islice)
483 #endif
484 cc        if (dyn_ss) then
485 cc        read(ientout,rec=i,err=101)
486 cc     &    ((csingle(l,k),l=1,3),k=1,nres),
487 cc     &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
488 cc     &    nss,(idssb(k),jdssb(k),k=1,nss),
489 cc     &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
490 cc        idssb(k)=idssb(k)-nres
491 cc        jdssb(k)=jdssb(k)-nres
492 cc        else
493         read(ientout,rec=i,err=101)
494      &    ((csingle(l,k),l=1,3),k=1,nres),
495      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
496      &    nss,(ihpb(k),jhpb(k),k=1,nss),
497      &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
498 cc         endif
499 cc        write (iout,*) 'CC', iR,ib,iparm,eini,efree
500         do j=1,2*nres
501           do k=1,3
502             c(k,j)=csingle(k,j)
503           enddo
504         enddo
505         call int_from_cart1(.false.)
506         iscore=0
507         if (indpdb.gt.0) then
508           call conf_compar(i,.false.,.true.)
509         endif 
510 c        if (dyn_ss) then
511         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)  
512      &    ((csingle(l,k),l=1,3),k=1,nres),
513      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
514      &    nss,(ihpb(k),jhpb(k),k=1,nss),
515 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore 
516      &    potE(i,nparmset),-entfac(i),rms_nat,iscore 
517 c        else
518           if (bxfile .or.cxfile .or. ensembles.gt.0) write
519      &     (ientin,rec=i)
520      &    ((csingle(l,k),l=1,3),k=1,nres),
521      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
522      &    nss,(ihpb(k),jhpb(k),k=1,nss),
523 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
524      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
525 c        endif
526 c        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
527 #ifndef MPI
528         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
529      &    -entfac(i),rms_nat,iscore)
530 #endif
531       enddo
532       close(ientout,status="delete")
533       close(istat)
534       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
535 #ifdef MPI
536       call MPI_Barrier(WHAM_COMM,IERROR)
537       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile 
538      &   .and. ensembles.eq.0) return
539       write (iout,*)
540       if (bxfile .or. ensembles.gt.0) then
541         if (nslice.eq.1) then
542           if (.not.separate_parset) then
543             bxname = prefix(:ilen(prefix))//".bx"
544           else
545             write (licz,'(bz,i3.3)') myparm
546             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
547           endif
548         else
549           if (.not.separate_parset) then
550             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
551           else
552             write (licz,'(bz,i3.3)') myparm
553             bxname = prefix(:ilen(prefix))//"par_"//licz//
554      &        "_slice_"//licz2//".bx"
555           endif
556         endif
557         open (ientout,file=bxname,status="unknown",
558      &      form="unformatted",access="direct",recl=lenrec1)
559         write (iout,*) "Master is creating binary database ",
560      &   bxname(:ilen(bxname))
561       endif
562       if (cxfile) then
563         if (nslice.eq.1) then
564           if (.not.separate_parset) then
565             cxname = prefix(:ilen(prefix))//".cx"
566           else
567             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
568           endif
569         else
570           if (.not.separate_parset) then
571             cxname = prefix(:ilen(prefix))//
572      &             "_slice_"//licz2//".cx"
573           else
574             cxname = prefix(:ilen(prefix))//"_par"//licz//
575      &             "_slice_"//licz2//".cx"
576           endif
577         endif
578 #if (defined(AIX) && !defined(JUBL))
579         call xdrfopen_(ixdrf,cxname, "w", iret)
580 #else
581         call xdrfopen(ixdrf,cxname, "w", iret)
582 #endif
583         if (iret.eq.0) then
584           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
585           cxfile=.false.
586         endif
587       endif
588       do j=0,nprocs-1
589         write (liczba,'(bz,i3.3)') j
590         if (separate_parset) then
591           write (licz,'(bz,i3.3)') myparm
592           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
593         else
594           bxname = prefix(:ilen(prefix))//liczba//".bx"
595         endif
596         open (ientin,file=bxname,status="unknown",
597      &    form="unformatted",access="direct",recl=lenrec1)
598         write (iout,*) "Master is reading conformations from ",
599      &   bxname(:ilen(bxname))
600         iii = 0
601 c        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
602 c        call flush(iout)
603         do i=indstart(j),indend(j)
604           iii = iii+1
605 cc          if (dyn_ss) then
606 cc          read(ientin,rec=iii,err=101)
607 cc     &      ((csingle(l,k),l=1,3),k=1,nres),
608 cc     &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
609 cc     &      nss,(idssb(k),jdssb(k),k=1,nss),
610 cc    &      eini,efree,rmsdev,iscor
611 cc          idssb(k)=idssb(k)-nres
612 cc          jdssb(k)=jdssb(k)-nres
613 cc          else
614           read(ientin,rec=iii,err=101)
615      &      ((csingle(l,k),l=1,3),k=1,nres),
616      &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
617      &      nss,(ihpb(k),jhpb(k),k=1,nss),
618      &      eini,efree,rmsdev,iscor
619 cc          endif
620           if (bxfile .or. ensembles.gt.0) then
621 cc          if (dyn_ss) then
622 cc            write (ientout,rec=i)
623 cc     &        ((csingle(l,k),l=1,3),k=1,nres),
624 cc     &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
625 cc     &        nss,(idssb(k)+nres,jdssb(k)+nres,k=1,nss),
626 cc    &        eini,efree,rmsdev,iscor
627 cc           else
628                         write (ientout,rec=i)
629      &        ((csingle(l,k),l=1,3),k=1,nres),
630      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
631      &        nss,(ihpb(k),jhpb(k),k=1,nss),
632      &        eini,efree,rmsdev,iscor
633 cc           write(iout,*) "W poszukiwaniu zlotych galotow"
634 cc           write(iout,*) "efree=",efree,iii
635 cc           endif
636           endif
637           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
638 #ifdef DEBUG
639           do k=1,2*nres
640             do l=1,3
641               c(l,k)=csingle(l,k)
642             enddo
643           enddo
644           call int_from_cart1(.false.)
645           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
646           write (iout,*) "The Cartesian geometry is:"
647           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
648           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
649           write (iout,*) "The internal geometry is:"
650           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
651           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
652           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
653           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
654           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
655           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
656           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
657           write (iout,'(f10.5,i5)') rmsdev,iscor
658 #endif
659         enddo ! i
660         write (iout,*) iii," conformations (from",indstart(j)," to",
661      &   indend(j),") read from ",
662      &   bxname(:ilen(bxname))
663         close (ientin,status="delete")
664       enddo ! j
665       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
666 #if (defined(AIX) && !defined(JUBL))
667       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
668 #else
669       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
670 #endif
671 #endif
672       return
673   101 write (iout,*) "Error in scratchfile."
674       call flush(iout)
675       return1
676       end
677 c-------------------------------------------------------------------------------
678       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
679       implicit none
680       include "DIMENSIONS"
681       include "DIMENSIONS.ZSCOPT"
682       include "DIMENSIONS.FREE"
683       include "DIMENSIONS.COMPAR"
684 #ifdef MPI
685       include "mpif.h"
686       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
687       include "COMMON.MPI"
688 #endif
689       include "COMMON.CONTROL"
690       include "COMMON.CHAIN"
691       include "COMMON.IOUNITS"
692       include "COMMON.PROTFILES"
693       include "COMMON.NAMES"
694       include "COMMON.VAR"
695       include "COMMON.SBRIDGE"
696       include "COMMON.GEO"
697       include "COMMON.FFIELD"
698       include "COMMON.ENEPS"
699       include "COMMON.LOCAL"
700       include "COMMON.WEIGHTS"
701       include "COMMON.INTERACT"
702       include "COMMON.FREE"
703       include "COMMON.ENERGIES"
704       include "COMMON.COMPAR"
705       include "COMMON.PROT"
706       integer i,j,itmp,iscor,iret,ixdrf
707       double precision rmsdev,efree,eini
708       real*4 csingle(3,maxres2),xoord(3,maxres2+2)
709       real*4 prec
710
711 c      write (iout,*) "cxwrite"
712 c      call flush(iout)
713       prec=10000.0
714       do i=1,nres
715        do j=1,3
716         xoord(j,i)=csingle(j,i)
717        enddo
718       enddo
719       do i=nnt,nct
720        do j=1,3
721         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
722        enddo
723       enddo
724
725       itmp=nres+nct-nnt+1
726
727 c      write (iout,*) "itmp",itmp
728 c      call flush(iout)
729 c       write (iout,*) "CNZ",eini,dyn_ss
730 #if (defined(AIX) && !defined(JUBL))
731       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
732
733 c      write (iout,*) "xdrf3dfcoord"
734 c      call flush(iout)
735       call xdrfint_(ixdrf, nss, iret)
736       do j=1,nss
737 cc       if (dyn_ss) then
738 cc       call xdrfint_(ixdrf, idssb(j)+nres, iret)
739 cc        call xdrfint_(ixdrf, jdssb(j)+nres, iret)
740 cc       else
741         call xdrfint_(ixdrf, ihpb(j), iret)
742         call xdrfint_(ixdrf, jhpb(j), iret)
743 cc       endif
744       enddo
745       call xdrffloat_(ixdrf,real(eini),iret) 
746       call xdrffloat_(ixdrf,real(efree),iret) 
747       call xdrffloat_(ixdrf,real(rmsdev),iret) 
748       call xdrfint_(ixdrf,iscor,iret) 
749 #else
750       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
751
752       call xdrfint(ixdrf, nss, iret)
753       do j=1,nss
754 cc       if (dyn_ss) then
755 cc        call xdrfint(ixdrf, idssb(j), iret)
756 cc        call xdrfint(ixdrf, jdssb(j), iret)
757 cc        idssb(j)=idssb(j)-nres
758 cc       jdssb(j)=jdssb(j)-nres
759 cc       else
760         call xdrfint(ixdrf, ihpb(j), iret)
761         call xdrfint(ixdrf, jhpb(j), iret)
762 cc       endif
763       enddo
764       call xdrffloat(ixdrf,real(eini),iret) 
765       call xdrffloat(ixdrf,real(efree),iret)
766       call xdrffloat(ixdrf,real(rmsdev),iret) 
767       call xdrfint(ixdrf,iscor,iret) 
768 #endif
769
770       return
771       end
772 c------------------------------------------------------------------------------
773       logical function conf_check(ii,iprint)
774       implicit none
775       include "DIMENSIONS"
776       include "DIMENSIONS.ZSCOPT"
777       include "DIMENSIONS.FREE"
778 #ifdef MPI
779       include "mpif.h"
780       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
781       include "COMMON.MPI"
782 #endif
783       include "COMMON.CHAIN"
784       include "COMMON.IOUNITS"
785       include "COMMON.PROTFILES"
786       include "COMMON.NAMES"
787       include "COMMON.VAR"
788       include "COMMON.SBRIDGE"
789       include "COMMON.GEO"
790       include "COMMON.FFIELD"
791       include "COMMON.ENEPS"
792       include "COMMON.LOCAL"
793       include "COMMON.WEIGHTS"
794       include "COMMON.INTERACT"
795       include "COMMON.FREE"
796       include "COMMON.ENERGIES"
797       include "COMMON.CONTROL"
798       include "COMMON.TORCNSTR"
799       integer j,k,l,ii,itj,iprint
800       if (.not.check_conf) then
801         conf_check=.true.
802         return
803       endif
804       call int_from_cart1(.false.)
805       do j=nnt+1,nct
806         if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
807           if (iprint.gt.0) 
808      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
809      &      " for conformation",ii
810           if (iprint.gt.1) then
811             write (iout,*) "The Cartesian geometry is:"
812             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
813             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
814             write (iout,*) "The internal geometry is:"
815             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
816             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
817             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
818             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
819             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
820             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
821           endif
822           if (iprint.gt.0) write (iout,*) 
823      &      "This conformation WILL NOT be added to the database."
824           conf_check=.false.
825           return
826         endif
827       enddo
828       do j=nnt,nct
829         itj=itype(j)
830         if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
831           if (iprint.gt.0) 
832      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
833      &     " for conformation",ii
834           if (iprint.gt.1) then
835             write (iout,*) "The Cartesian geometry is:"
836             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
837             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
838             write (iout,*) "The internal geometry is:"
839             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
840             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
841             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
842             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
843             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
844             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
845           endif
846           if (iprint.gt.0) write (iout,*) 
847      &      "This conformation WILL NOT be added to the database."
848           conf_check=.false.
849           return
850         endif
851       enddo
852       do j=3,nres
853         if (theta(j).le.0.0d0) then
854           if (iprint.gt.0) 
855      &    write (iout,*) "Zero theta angle(s) in conformation",ii
856           if (iprint.gt.1) then
857             write (iout,*) "The Cartesian geometry is:"
858             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
859             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
860             write (iout,*) "The internal geometry is:"
861             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
862             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
863             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
864             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
865             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
866             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
867           endif
868           if (iprint.gt.0) write (iout,*)
869      &      "This conformation WILL NOT be added to the database." 
870           conf_check=.false.
871           return
872         endif
873         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
874       enddo
875       conf_check=.true.
876 c      write (iout,*) "conf_check passed",ii
877       return
878       end