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