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