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