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