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