source/wham/io_database.f90

   1       module io_database
   2 !-----------------------------------------------------------------------------
   3       use names
   4       use wham_data
   5       use io_units
   6       use io_base, only:ilen
   7       use energy_data, only:nnt,nct,nss,ihpb,jhpb
   8       use MD_data, only:iset
   9       use geometry_data, only:nres,c
  10 #ifdef MPI
  11       use MPI_data
  12 !      include "COMMON.MPI"
  13 #endif
  14
  15       implicit none
  16 !-----------------------------------------------------------------------------
  17 !
  18 !
  19 !-----------------------------------------------------------------------------
  20       contains
  21 !-----------------------------------------------------------------------------
  22 ! readrtns.F
  23 !-------------------------------------------------------------------------------
  24       subroutine opentmp(islice,iunit,bprotfile_temp)
  25 !      implicit none
  26 !      include "DIMENSIONS"
  27 !      include "DIMENSIONS.ZSCOPT"
  28 !      include "DIMENSIONS.FREE"
  29 !      use MPI_data, only:me
  30 #ifdef MPI
  31       include "mpif.h"
  32       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
  33 !      include "COMMON.MPI"
  34 #endif
  35 !      include "COMMON.IOUNITS"
  36 !      include "COMMON.PROTFILES"
  37 !      include "COMMON.PROT"
  38 !      include "COMMON.FREE"
  39       character(len=64) :: bprotfile_temp
  40       character(len=3) :: liczba,liczba2
  41       character(len=2) :: liczba1
  42       integer :: iunit,islice
  43 !      integer ilen,iroof
  44 !      external ilen,iroof
  45 !      logical :: lerr
  46 !      integer :: lenrec,lenrec2
  47
  48 !el
  49 !      lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
  50 !      lenrec=lenrec2+8
  51       write (liczba1,'(bz,i2.2)') islice
  52 #ifdef MPI
  53       write (liczba,'(bz,i3.3)') me
  54 !#ifdef MPI
  55 !      write (iout,*) "separate_parset ",separate_parset,
  56 !     &  " myparm",myparm
  57       if (separate_parset) then
  58       write (liczba2,'(bz,i3.3)') myparm
  59       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"// &
  60         prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
  61       open (iunit,file=bprotfile_temp,status="unknown",&
  62           form="unformatted",access="direct",recl=lenrec)
  63       else
  64       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"// &
  65         prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
  66       open (iunit,file=bprotfile_temp,status="unknown",&
  67           form="unformatted",access="direct",recl=lenrec)
  68       endif
  69 #else
  70       bprotfile_temp = scratchdir(:ilen(scratchdir))// &
  71         "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
  72       open (iunit,file=bprotfile_temp,status="unknown",&
  73           form="unformatted",access="direct",recl=lenrec)
  74 #endif
  75 !      write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
  76 !     &  bprotfile_temp
  77 !      call flush(iout)
  78       return
  79       end subroutine opentmp
  80 !-------------------------------------------------------------------------------
  81       subroutine read_database(*)
  82
  83 !      use energy_data, only:nct,nnt,nss
  84 !      implicit none
  85 !      include "DIMENSIONS"
  86 !      include "DIMENSIONS.ZSCOPT"
  87 !      include "DIMENSIONS.FREE"
  88       use MPI_data, only:me,nprocs
  89 #ifdef MPI
  90       include "mpif.h"
  91       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
  92 !      include "COMMON.MPI"
  93 #endif
  94 !      include "COMMON.CHAIN"
  95 !      include "COMMON.IOUNITS"
  96 !      include "COMMON.PROTFILES"
  97 !      include "COMMON.NAMES"
  98 !      include "COMMON.VAR"
  99 !      include "COMMON.GEO"
 100 !      include "COMMON.ENEPS"
 101 !      include "COMMON.PROT"
 102 !      include "COMMON.INTERACT"
 103 !      include "COMMON.FREE"
 104 !      include "COMMON.SBRIDGE"
 105 !      include "COMMON.OBCINKA"
 106       real(kind=4) :: csingle(3,nres*2) !(3,maxres2)
 107       character(len=64) :: nazwa,bprotfile_temp
 108       character(len=3) :: liczba
 109       character(len=2) :: liczba1
 110       integer :: i,j,ii,jj(nslice),k,kk(nslice),l,&
 111         ll(nslice),mm(nslice),if
 112       integer :: nrec,nlines,iscor,iunit,islice
 113       real(kind=8) :: energ
 114 !      integer ilen,iroof
 115 !      external ilen,iroof
 116       real(kind=8) :: rmsdev,energia(0:n_ene),efree,eini,temp
 117 !el      real(kind=8) :: rmsdev,energia(0:max_eneW),efree,eini,temp
 118       real(kind=8) :: prop(nQ) !(maxQ)
 119       integer :: ntot_all(nslice,0:nprocs-1)!(maxslice,0:maxprocs-1)
 120       integer :: iparm,ib,iib,ir,nprop,nthr,npars
 121       real(kind=8) :: etot,time
 122       integer :: ixdrf,iret
 123       logical :: lerr,linit
 124
 125       lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
 126       lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
 127       lenrec=lenrec2+8
 128       write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,&
 129         " lenrec2",lenrec2
 130
 131       do i=1,nQ
 132         prop(i)=0.0d0
 133       enddo
 134       do islice=1,nslice
 135         ll(islice)=0
 136         mm(islice)=0
 137       enddo
 138       write (iout,*) "nparmset",nparmset
 139       if (hamil_rep) then
 140         npars=1
 141       else
 142         npars=nparmset
 143       endif
 144       do iparm=1,npars
 145
 146       if (replica(iparm)) then
 147         nthr = 1
 148       else
 149         nthr = nT_h(iparm)
 150       endif
 151
 152       do ib=1,nthr
 153       do iR=1,nRR(ib,iparm)
 154
 155       write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 156       do islice=1,nslice
 157         jj(islice)=0
 158         kk(islice)=0
 159       enddo
 160
 161       IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
 162 ! Read conformations from binary DA files (one per batch) and write them to
 163 ! a binary DA scratchfile.
 164         write (liczba,'(bz,i3.3)') me
 165         do if=1,nfile_bin(iR,ib,iparm)
 166           nazwa=protfiles(if,1,iR,ib,iparm) &
 167            (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
 168           open (ientin,file=nazwa,status="old",form="unformatted",&
 169            access="direct",recl=lenrec2,err=1111)
 170           ii=0
 171           do islice=1,nslice
 172             call opentmp(islice,ientout,bprotfile_temp)
 173             call bxread(nazwa,islice,ii,jj(islice),kk(islice),ll(islice),&
 174               mm(islice),iR,ib,iparm)
 175             close(ientout)
 176           enddo
 177           close(ientin)
 178         enddo
 179       ENDIF ! NFILE_BIN>0
 180 !
 181       IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
 182 ! Read conformations from multiple ASCII int files and write them to a binary
 183 ! DA scratchfile.
 184         do if=1,nfile_asc(iR,ib,iparm)
 185           nazwa=protfiles(if,2,iR,ib,iparm) &
 186            (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
 187           open(unit=ientin,file=nazwa,status='old',err=1111)
 188           write(iout,*) "reading ",nazwa(:ilen(nazwa))
 189           ii=0
 190           call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
 191         enddo ! if
 192       ENDIF
 193       IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
 194 ! Read conformations from cx files and write them to a binary
 195 ! DA scratchfile.
 196         do if=1,nfile_cx(iR,ib,iparm)
 197           nazwa=protfiles(if,2,iR,ib,iparm) &
 198            (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
 199           write(iout,*) "reading ",nazwa(:ilen(nazwa))
 200           ii=0
 201           print *,"Calling cxread"
 202           call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,&
 203              *1111)
 204 write(iout,*)"after call cxread"
 205           close(ientout)
 206           write (iout,*) "exit cxread"
 207           call flush(iout)
 208         enddo
 209       ENDIF
 210 write(iout,*)"*********************in read database"
 211
 212       do islice=1,nslice
 213 !        stot(islice)=0
 214         stot(islice)=stot(islice)+jj(islice)
 215       enddo
 216
 217       enddo
 218       enddo
 219       write (iout,*) "IPARM",iparm
 220       enddo
 221
 222       if (nslice.eq.1) then
 223 #ifdef MPI
 224         write (liczba,'(bz,i3.3)') me
 225         bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"// &
 226           prefix(:ilen(prefix))//liczba//".xbin.tmp"
 227 #else
 228         bprotfile_temp = scratchdir(:ilen(scratchdir))// &
 229            "/"//prefix(:ilen(prefix))//".xbin.tmp"
 230 #endif
 231         write(iout,*) mm(1)," conformations read",ll(1),&
 232           " conformations written to ",&
 233           bprotfile_temp(:ilen(bprotfile_temp))
 234       else
 235         do islice=1,nslice
 236           write (liczba1,'(bz,i2.2)') islice
 237 #ifdef MPI
 238           write (liczba,'(bz,i3.3)') me
 239           bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"// &
 240             prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
 241 #else
 242           bprotfile_temp = scratchdir(:ilen(scratchdir))// &
 243              "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
 244 #endif
 245           write(iout,*) mm(islice)," conformations read",ll(islice),&
 246           " conformations written to ",&
 247           bprotfile_temp(:ilen(bprotfile_temp))
 248         enddo
 249       endif
 250
 251 #ifdef MPI
 252 ! Check if everyone has the same number of conformations
 253       call MPI_Allgather(stot(1),nslice,MPI_INTEGER,&
 254         ntot_all(1,0),nslice,MPI_INTEGER,MPI_Comm_World,IERROR)
 255       lerr=.false.
 256       do i=0,nprocs-1
 257         if (i.ne.me) then
 258           do islice=1,nslice
 259           if (stot(islice).ne.ntot_all(islice,i)) then
 260             write (iout,*) "Number of conformations at processor",i,&
 261              " differs from that at processor",me,&
 262              stot(islice),ntot_all(islice,i)," slice",islice
 263             lerr = .true.
 264           endif
 265           enddo
 266         endif
 267       enddo
 268       if (lerr) then
 269         write (iout,*)
 270         write (iout,*) "Numbers of conformations read by processors"
 271         write (iout,*)
 272         do i=0,nprocs-1
 273           write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
 274         enddo
 275         write (iout,*) "Calculation terminated."
 276         call flush(iout)
 277         return 1
 278       endif
 279       do islice=1,nslice
 280         ntot(islice)=stot(islice)
 281       enddo
 282 write(iout,*) "end of read database"
 283       return
 284 #endif
 285  1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
 286       call flush(iout)
 287       return 1
 288       end subroutine read_database
 289 !--------------------------------------------------------------------------------
 290       integer function iroof(n,m)
 291       integer :: n,m,ii
 292       ii = n/m
 293       if (ii*m .lt. n) ii=ii+1
 294       iroof = ii
 295       return
 296       end function iroof
 297 !--------------------------------------------------------------------------------
 298 ! bxread.F
 299 !--------------------------------------------------------------------------------
 300       subroutine bxread(nazwa,islice,ii,jj,kk,ll,mm,iR,ib,iparm)
 301 !      implicit none
 302 !      include "DIMENSIONS"
 303 !      include "DIMENSIONS.ZSCOPT"
 304 !      include "DIMENSIONS.FREE"
 305 !      use energy_data, only:nnt,nct,nss,ihpb,jhpbi
 306       use MPI_data, only:nprocs
 307 #ifdef MPI
 308       include "mpif.h"
 309       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 310 !      include "COMMON.MPI"
 311 #endif
 312 !      include "COMMON.CHAIN"
 313 !      include "COMMON.IOUNITS"
 314 !      include "COMMON.PROTFILES"
 315 !      include "COMMON.NAMES"
 316 !      include "COMMON.VAR"
 317 !      include "COMMON.GEO"
 318 !      include "COMMON.ENEPS"
 319 !      include "COMMON.PROT"
 320 !      include "COMMON.INTERACT"
 321 !      include "COMMON.FREE"
 322 !      include "COMMON.SBRIDGE"
 323       real(kind=4) :: csingle(3,nres*2) !(3,maxres2)
 324       character(len=64) :: nazwa,bprotfile_temp
 325       character(len=3) :: liczba
 326       integer :: i,is,ie,j,ii,jj,k,kk,l,ll,mm,if
 327       integer :: nrec,nlines,iscor,islice
 328       real(kind=8) :: energ
 329 !      integer ilen,iroof
 330 !      external ilen,iroof
 331       real(kind=8) :: rmsdev,energia(0:n_ene),efree,eini,temp
 332 !el      real(kind=8) :: rmsdev,energia(0:max_eneW),efree,eini,temp
 333       real(kind=8) :: prop(nQ) !(maxQ)
 334       integer :: ntot_all(0:nprocs-1)!(0:maxprocs-1)
 335       integer :: iparm,ib,iib,ir,nprop,nthr,nrec_slice
 336       real(kind=8) :: etot,time
 337       logical :: lerr
 338       nrec_slice=(rec_end(iR,ib,iparm)-rec_start(iR,ib,iparm)+1)/nslice
 339       is=rec_start(iR,ib,iparm)+(islice-1)*nrec_slice
 340       ie=rec_start(iR,ib,iparm)+islice*nrec_slice-1
 341       write (iout,*) "bxread: islice",islice," nslice",nslice,&
 342        " nrec_slice",nrec_slice
 343       write (iout,*) "is",is," ie",ie,"rec_start",&
 344         rec_start(iR,ib,iparm)," rec_end",rec_end(iR,ib,iparm)
 345       do i=is,ie
 346             read(ientin,rec=i+1,err=101) &
 347               ((csingle(l,k),l=1,3),k=1,nres),&
 348               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 349               nss,(ihpb(k),jhpb(k),k=1,nss),&
 350               eini,efree,rmsdev,(prop(j),j=1,nQ),iscor
 351             ii=ii+1
 352             kk=kk+1
 353             if (mod(kk,isampl(iparm)).eq.0) then
 354             jj=jj+1
 355             write(ientout,rec=jj) &
 356               ((csingle(l,k),l=1,3),k=1,nres),&
 357               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 358               nss,(ihpb(k),jhpb(k),k=1,nss),&
 359               eini,efree,rmsdev,(prop(j),j=1,nQ),iR,ib,iparm
 360 #ifdef DEBUG
 361             do i=1,2*nres
 362               do j=1,3
 363                 c(j,i)=csingle(j,i)
 364               enddo
 365             enddo
 366             call int_from_cart1(.false.)
 367             write (iout,*) "Writing conformation, record",jj
 368             write (iout,*) "Cartesian coordinates"
 369             write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 370             write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 371             write (iout,*) "Internal coordinates"
 372             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 373             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 374             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 375             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 376             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 377             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 378             write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 379             write (iout,'(f10.5,i5)') rmsdev,iscor
 380 #endif
 381             endif
 382           enddo
 383   101     continue
 384           close(ientin)
 385           write (iout,*) ii," conformations read from DA file ",&
 386             nazwa(:ilen(nazwa))
 387           write (iout,*) kk," conformations read so far, slice",islice
 388           write (iout,*) jj," conformations stored so far, slice",islice
 389
 390       return
 391       end subroutine bxread
 392 !--------------------------------------------------------------------------------
 393 ! cxread.F
 394 !--------------------------------------------------------------------------------
 395       subroutine cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,*)
 396
 397 #define DEBUG
 398 #ifdef DEBUG
 399       use geometry, only:int_from_cart1
 400       use geometry_data, only:vbld,rad2deg,theta,phi,alph,omeg
 401       integer :: iscor
 402 #endif
 403 !      implicit real*8 (a-h,o-z)
 404 !      include 'DIMENSIONS'
 405 !      include 'DIMENSIONS.ZSCOPT'
 406 !      include 'DIMENSIONS.FREE'
 407       integer,parameter :: MaxTraj=2050
 408 !      include 'COMMON.CHAIN'
 409 !      include 'COMMON.INTERACT'
 410 !      include 'COMMON.NAMES'
 411 !      include 'COMMON.IOUNITS'
 412 !      include 'COMMON.HEADER'
 413 !      include 'COMMON.SBRIDGE'
 414 !      include 'COMMON.PROTFILES'
 415 !      include 'COMMON.OBCINKA'
 416 !      include 'COMMON.FREE'
 417 !      include 'COMMON.VAR'
 418 !      include 'COMMON.GEO'
 419 !      include 'COMMON.PROT'
 420       character(len=64) :: nazwa,bprotfile_temp
 421       real(kind=4) :: rtime,rpotE,ruconst,rt_bath,rprop(nQ) !(2000) !(maxQ)
 422       real(kind=8) :: time
 423       integer :: iret,itmp,itraj,ntraj
 424       real(kind=4) :: xoord(3,2*nres+2),prec
 425       integer :: nstep(0:MaxTraj-1)
 426 !      integer ilen
 427 !      external ilen
 428       integer :: ii,jj(nslice),kk(nslice),ll(nslice),mm(nslice) !(maxslice)
 429       integer :: is(nSlice),ie(nSlice),nrec_slice
 430       real(kind=8) :: ts(nSlice),te(nSlice),time_slice
 431       integer :: iR,ib,iparm,i,j,it,islice,nprop_prev
 432       integer :: k,l,iib,islice1,nprop
 433       real(kind=8) :: efree,rmsdev
 434       integer :: ixdrf
 435 !el      integer :: slice
 436 !      logical :: conf_check
 437 !      ixdrf=0
 438 !      nprop=0
 439
 440 !      ruconst=0.0d0
 441 !      rtime=0.0d0
 442 !      rpotE=0.0d0
 443 !      rt_bath=0.0d0
 444
 445       call set_slices(is,ie,ts,te,iR,ib,iparm)
 446       nprop_prev=0
 447       do i=1,nQ
 448         rprop(i)=0.0d0
 449       enddo
 450       do i=0,MaxTraj-1
 451         nstep(i)=0
 452       enddo
 453       ntraj=0
 454       it=0
 455       iret=1
 456 #if (defined(AIX) && !defined(JUBL))
 457       call xdrfopen_(ixdrf,nazwa, "r", iret)
 458 #else
 459       call xdrfopen(ixdrf,nazwa, "r", iret)
 460 #endif
 461       if (iret.eq.0) return 1
 462
 463       islice1=1
 464       call opentmp(islice1,ientout,bprotfile_temp)
 465       print *,"bumbum" !d
 466       do while (iret.gt.0)
 467
 468 #if (defined(AIX) && !defined(JUBL))
 469       call xdrffloat_(ixdrf, rtime, iret)
 470       print *,"rtime",rtime," iret",iret !d
 471       call xdrffloat_(ixdrf, rpotE, iret)
 472       write (iout,*) "rpotE",rpotE," iret",iret !d
 473       call flush(iout)
 474       call xdrffloat_(ixdrf, ruconst, iret)
 475       call xdrffloat_(ixdrf, rt_bath, iret)
 476       call xdrfint_(ixdrf, nss, iret)
 477       do j=1,nss
 478         call xdrfint_(ixdrf, ihpb(j), iret)
 479         call xdrfint_(ixdrf, jhpb(j), iret)
 480       enddo
 481       call xdrfint_(ixdrf, nprop, iret)
 482       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep) &
 483         call xdrfint(ixdrf, iset, iret)
 484       do i=1,nprop
 485         call xdrffloat_(ixdrf, rprop(i), iret)
 486       enddo
 487 #else
 488       call xdrffloat(ixdrf, rtime, iret)
 489       call xdrffloat(ixdrf, rpotE, iret)
 490       write (iout,*) "rpotE",rpotE," iret",iret !d
 491       call flush(iout)
 492       call xdrffloat(ixdrf, ruconst, iret)
 493       call xdrffloat(ixdrf, rt_bath, iret)
 494       call xdrfint(ixdrf, nss, iret)
 495       do j=1,nss
 496         call xdrfint(ixdrf, ihpb(j), iret)
 497         call xdrfint(ixdrf, jhpb(j), iret)
 498       enddo
 499       call xdrfint(ixdrf, nprop, iret)
 500       write (iout,*) "nprop",nprop !d
 501       if (it.gt.0 .and. nprop.ne.nprop_prev) then
 502         write (iout,*) "Warning previous nprop",nprop_prev,&
 503          " current",nprop
 504         nprop=nprop_prev
 505       else
 506         nprop_prev=nprop
 507       endif
 508       call flush(iout)
 509       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep) &
 510         call xdrfint(ixdrf, iset, iret)
 511       do i=1,nprop
 512         call xdrffloat(ixdrf, rprop(i), iret)
 513       enddo
 514 #endif
 515       if (iret.eq.0) exit
 516       itraj=mod(it,totraj(iR,iparm))
 517 #define DEBUG
 518 #ifdef DEBUG
 519       write (iout,*) "ii",ii," itraj",itraj," it",it
 520 #endif
 521       if (iset.eq.0) iset = 1
 522       call flush(iout)
 523       it=it+1
 524       if (itraj.gt.ntraj) ntraj=itraj
 525       nstep(itraj)=nstep(itraj)+1
 526 !      rprop(2)=dsqrt(rprop(2))
 527 !      rprop(3)=dsqrt(rprop(3))
 528 #ifdef DEBUG
 529        write (iout,*) "umbrella ",umbrella
 530        write (iout,*) rtime,rpotE,rt_bath,nss,&
 531            (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
 532        write (iout,*) "nprop",nprop," iset",iset," myparm",myparm
 533        call flush(iout)
 534 #endif
 535       prec=10000.0
 536       itmp=0
 537 #if (defined(AIX) && !defined(JUBL))
 538       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 539 #else
 540       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 541 #endif
 542 #ifdef DEBUG
 543       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,2*nres+2)
 544 #endif
 545 #undef DEBUG
 546       if (iret.eq.0) exit
 547       if (itmp .ne. nres + nct - nnt + 1) then
 548         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
 549         call flush(iout)
 550         exit
 551       endif
 552
 553       time=rtime
 554       write (iout,*) "calling slice" !d
 555       call flush(iout) !d
 556       islice=slice(nstep(itraj),time,is,ie,ts,te)
 557       write (iout,*) "islice",islice !d
 558       call flush(iout) !d
 559
 560       do i=1,nres
 561         do j=1,3
 562           c(j,i)=xoord(j,i)
 563         enddo
 564       enddo
 565       do i=1,nct-nnt+1
 566         do j=1,3
 567           c(j,i+nres+nnt-1)=xoord(j,i+nres)
 568         enddo
 569       enddo
 570
 571       if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset &
 572           .or. iset.eq.myparm)) then
 573         ii=ii+1
 574         kk(islice)=kk(islice)+1
 575         mm(islice)=mm(islice)+1
 576         if (mod(nstep(itraj),isampl(iparm)).eq.0 .and. &
 577            conf_check(ll(islice)+1,1)) then
 578           if (replica(iparm)) then
 579              rt_bath=1.0d0/(rt_bath*1.987D-3)
 580              do i=1,nT_h(iparm)
 581                if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
 582                  iib = i
 583                  goto 22
 584                endif
 585              enddo
 586   22         continue
 587              if (i.gt.nT_h(iparm)) then
 588                write (iout,*) "Error - temperature of conformation",&
 589                ii,1.0d0/(rt_bath*1.987D-3),&
 590                " does not match any of the list"
 591                write (iout,*) &
 592                 1.0d0/(rt_bath*1.987D-3),&
 593                 (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 594                call flush(iout)
 595 !               exit
 596 !               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 597                ii=ii-1
 598                kk(islice)=kk(islice)-1
 599                mm(islice)=mm(islice)-1
 600                goto 112
 601              endif
 602           else
 603             iib = ib
 604           endif
 605
 606           efree=0.0d0
 607           jj(islice)=jj(islice)+1
 608           if (umbrella(iparm)) then
 609             snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 610           else if (hamil_rep) then
 611             snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
 612           else
 613             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 614           endif
 615           ll(islice)=ll(islice)+1
 616 #ifdef DEBUG
 617           write (iout,*) "Writing conformation, record",ll(islice)
 618           write (iout,*) "ib",ib," iib",iib
 619           write (iout,*) "ntraj",ntraj," itraj",itraj,&
 620             " nstep",nstep(itraj)
 621           write (iout,*) "pote",rpotE," time",rtime
 622 !          if (replica(iparm)) then
 623 !            write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
 624 !            write (iout,*) "TEMP list"
 625 !            write (iout,*)
 626 !     &       (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 627 !          endif
 628           write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 629 !          write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 630 !          write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 631           call flush(iout)
 632 #endif
 633           if (islice.ne.islice1) then
 634 !            write (iout,*) "islice",islice," islice1",islice1
 635             close(ientout)
 636 !            write (iout,*) "Closing file ",
 637 !     &          bprotfile_temp(:ilen(bprotfile_temp))
 638             call opentmp(islice,ientout,bprotfile_temp)
 639 !            write (iout,*) "Opening file ",
 640 !     &          bprotfile_temp(:ilen(bprotfile_temp))
 641             islice1=islice
 642           endif
 643           if (umbrella(iparm)) then
 644             write(ientout,rec=ll(islice)) &
 645               ((xoord(l,k),l=1,3),k=1,nres),&
 646               ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),&
 647               nss,(ihpb(k),jhpb(k),k=1,nss),&
 648               rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),&
 649               iset,iib,iparm
 650           else if (hamil_rep) then
 651             write(ientout,rec=ll(islice)) &
 652               ((xoord(l,k),l=1,3),k=1,nres),&
 653               ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),&
 654               nss,(ihpb(k),jhpb(k),k=1,nss),&
 655               rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),&
 656               iR,iib,iset
 657           else
 658             write(ientout,rec=ll(islice)) &
 659               ((xoord(l,k),l=1,3),k=1,nres),&
 660               ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),&
 661               nss,(ihpb(k),jhpb(k),k=1,nss),&
 662               rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),&
 663               iR,iib,iparm
 664           endif
 665 #ifdef DEBUG
 666           call int_from_cart1(.false.)
 667           write (iout,*) "Writing conformation, record",ll(islice)
 668           write (iout,*) "Cartesian coordinates"
 669           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 670           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 671           write (iout,*) "Internal coordinates"
 672           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 673           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 674           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 675           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 676           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 677           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 678           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 679 !          write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
 680           write (iout,'(16i5)') iscor
 681           call flush(iout)
 682 #endif
 683         endif
 684       endif
 685
 686   112 continue
 687
 688       enddo
 689       close(ientout)
 690 #if (defined(AIX) && !defined(JUBL))
 691       call xdrfclose_(ixdrf, iret)
 692 #else
 693       call xdrfclose(ixdrf, iret)
 694 #endif
 695       write (iout,'(i10," trajectories found in file.")') ntraj+1
 696       write (iout,'(a)') "Numbers of steps in trajectories:"
 697       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 698       write (iout,*) ii," conformations read from file",&
 699          nazwa(:ilen(nazwa))
 700       do islice=1,nslice
 701         write (iout,*) mm(islice)," conformations read so far, slice",&
 702           islice
 703         write (iout,*) ll(islice),&
 704         " conformations stored so far, slice",islice
 705       enddo
 706       call flush(iout)
 707 #undef DEBUG
 708       return
 709       end subroutine cxread
 710 !--------------------------------------------------------------------------------
 711 ! xread.F
 712 !--------------------------------------------------------------------------------
 713       subroutine xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
 714
 715       use geometry_data
 716 !      implicit none
 717 !      include "DIMENSIONS"
 718 !      include "DIMENSIONS.ZSCOPT"
 719 !      include "DIMENSIONS.FREE"
 720       use MPI_data, only:nprocs
 721 #ifdef MPI
 722       include "mpif.h"
 723       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 724 !      include "COMMON.MPI"
 725 #endif
 726       integer,parameter :: MaxTraj=2050
 727 !      include "COMMON.CHAIN"
 728 !      include "COMMON.IOUNITS"
 729 !      include "COMMON.PROTFILES"
 730 !      include "COMMON.NAMES"
 731 !      include "COMMON.VAR"
 732 !      include "COMMON.GEO"
 733 !      include "COMMON.ENEPS"
 734 !      include "COMMON.PROT"
 735 !      include "COMMON.INTERACT"
 736 !      include "COMMON.FREE"
 737 !      include "COMMON.SBRIDGE"
 738 !      include "COMMON.OBCINKA"
 739       real(kind=4) :: csingle(3,nres*2)
 740       character(len=64) :: nazwa,bprotfile_temp
 741       integer :: i,j,k,l,ii,jj(nslice),kk(nslice),ll(nslice),&
 742         mm(nslice) !(maxslice)
 743       integer :: iscor,islice,islice1 !el,slice
 744       real(kind=8) :: energ
 745 !      integer ilen,iroof
 746 !      external ilen,iroof
 747       real(kind=8) :: rmsdev,energia(0:n_ene),efree,eini,temp
 748 !el      real(kind=8) :: rmsdev,energia(0:max_eneW),efree,eini,temp
 749       real(kind=8) :: prop(nQ) !(maxQ)
 750       integer :: ntot_all(0:nprocs-1)!(0:maxprocs-1)
 751       integer :: iparm,ib,iib,ir,nprop,nthr
 752       real(kind=8) :: etot,time,ts(nslice),te(nslice)
 753       integer :: is(nslice),ie(nslice),itraj,ntraj,it,iset
 754       integer :: nstep(0:MaxTraj-1)
 755       logical :: lerr
 756
 757       call set_slices(is,ie,ts,te,iR,ib,iparm)
 758       do i=1,nQ
 759         prop(i)=0.0d0
 760       enddo
 761       do i=0,MaxTraj-1
 762         nstep(i)=0
 763       enddo
 764       ntraj=0
 765       it=0
 766       islice1=1
 767       call opentmp(islice1,ientout,bprotfile_temp)
 768       do while (.true.)
 769         if (replica(iparm)) then
 770           if (hamil_rep .or. umbrella(iparm)) then
 771           read (ientin,*,end=1112,err=1112) time,eini,&
 772             etot,temp,nss,(ihpb(j),jhpb(j),j=1,nss),&
 773             nprop,(prop(j),j=1,nprop),iset
 774           else
 775           read (ientin,*,end=1112,err=1112) time,eini,&
 776             etot,temp,nss,(ihpb(j),jhpb(j),j=1,nss),&
 777             nprop,(prop(j),j=1,nprop)
 778           endif
 779           temp=1.0d0/(temp*1.987D-3)
 780 !           write (iout,*) time,eini,etot,nss,
 781 !     &     (ihpb(j),jhpb(j),j=1,nss),(prop(j),j=1,nprop)
 782 !           call flush(iout)
 783            do i=1,nT_h(iparm)
 784              if (beta_h(i,iparm).eq.temp) then
 785                iib = i
 786                goto 22
 787              endif
 788            enddo
 789   22       continue
 790            if (i.gt.nT_h(iparm)) then
 791              write (iout,*) "Error - temperature of conformation",&
 792              ii,1.0d0/(temp*1.987D-3),&
 793              " does not match any of the list"
 794              write (iout,*) &
 795               1.0d0/(temp*1.987D-3),&
 796               (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 797              call flush(iout)
 798 #ifdef MPI
 799              call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 800 #endif
 801            endif
 802         else
 803            read (ientin,*,end=1112,err=1112) time,eini,&
 804              etot,nss,(ihpb(j),jhpb(j),j=1,nss),&
 805              nprop,(prop(j),j=1,nprop)
 806              iib = ib
 807         endif
 808         itraj=mod(it,totraj(iR,iparm))
 809 !        write (*,*) "ii",ii," itraj",itraj
 810 !        call flush(iout)
 811         it=it+1
 812         if (itraj.gt.ntraj) ntraj=itraj
 813         nstep(itraj)=nstep(itraj)+1
 814         islice=slice(nstep(itraj),time,is,ie,ts,te)
 815         read (ientin,'(8f10.5)',end=1112,err=1112) &
 816           ((csingle(l,k),l=1,3),k=1,nres),&
 817           ((csingle(l,k+nres),l=1,3),k=nnt,nct)
 818         efree=0.0d0
 819         if (islice.gt.0 .and. islice.le.nslice) then
 820         ii=ii+1
 821         kk(islice)=kk(islice)+1
 822         mm(islice)=mm(islice)+1
 823         if (mod(nstep(itraj),isampl(iparm)).eq.0) then
 824         jj(islice)=jj(islice)+1
 825         if (hamil_rep) then
 826           snk(iR,iib,iset,islice)=snk(iR,iib,iset,islice)+1
 827         else if (umbrella(iparm)) then
 828           snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 829         else
 830           snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 831         endif
 832         ll(islice)=ll(islice)+1
 833 !         write (iout,*) ii,kk,jj,ll,eini,(prop(j),j=1,nprop)
 834 #ifdef DEBUG
 835 !        write (iout,*) "Writing conformation, record",ll(islice)
 836 !        write (iout,*) "ib",ib," iib",iib
 837          if (replica(iparm)) then
 838            write (iout,*) "TEMP",1.0d0/(temp*1.987D-3)
 839            write (iout,*) "TEMP list"
 840            write (iout,*) &
 841             (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 842          endif
 843          call flush(iout)
 844 #endif
 845 !         write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 846 !         write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 847 !         write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 848 !         call flush(iout)
 849          if (islice.ne.islice1) then
 850 !             write (iout,*) "islice",islice," islice1",islice1
 851              close(ientout)
 852 !             write (iout,*) "Closing file ",
 853 !     &             bprotfile_temp(:ilen(bprotfile_temp))
 854              call opentmp(islice,ientout,bprotfile_temp)
 855 !             write (iout,*) "Opening file ",
 856 !     &             bprotfile_temp(:ilen(bprotfile_temp))
 857 !             call flush(iout)
 858              islice1=islice
 859          endif
 860          write(ientout,rec=ll(islice)) &
 861            ((csingle(l,k),l=1,3),k=1,nres),&
 862            ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 863            nss,(ihpb(k),jhpb(k),k=1,nss),&
 864            eini,efree,rmsdev,(prop(i),i=1,nQ),iR,iib,iparm
 865 #ifdef DEBUG
 866          do i=1,2*nres
 867            do j=1,3
 868              c(j,i)=csingle(j,i)
 869            enddo
 870          enddo
 871          call int_from_cart1(.false.)
 872          write (iout,*) "Writing conformation, record",ll(islice)
 873          write (iout,*) "Cartesian coordinates"
 874          write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 875          write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 876          write (iout,*) "Internal coordinates"
 877          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 878          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 879          write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 880          write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 881          write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 882          write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 883          write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 884 !         write (iout,'(8f10.5)') (prop(j),j=1,nQ)
 885          write (iout,'(16i5)') iscor
 886          call flush(iout)
 887 #endif
 888          endif
 889          endif
 890        enddo
 891  1112  continue
 892        close(ientout)
 893        write (iout,'(i10," trajectories found in file.")') ntraj+1
 894        write (iout,'(a)') "Numbers of steps in trajectories:"
 895        write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 896        write (iout,*) ii," conformations read from file",&
 897          nazwa(:ilen(nazwa))
 898        write (iout,*) mm(islice)," conformations read so far, slice",&
 899           islice
 900        write (iout,*) ll(islice)," conformations stored so far, slice",&
 901          islice
 902        call flush(iout)
 903        return
 904        end subroutine xread
 905 !--------------------------------------------------------------------------------
 906 ! enecalc1.F
 907 !--------------------------------------------------------------------------------
 908       subroutine write_dbase(islice,*)
 909
 910       use geometry_data
 911       use control_data, only:indpdb
 912       use w_compar_data
 913       use conform_compar, only:conf_compar
 914 !      implicit none
 915 !      include "DIMENSIONS"
 916 !      include "DIMENSIONS.ZSCOPT"
 917 !      include "DIMENSIONS.FREE"
 918 !      include "DIMENSIONS.COMPAR"
 919       use geometry, only:int_from_cart1
 920 #ifdef MPI
 921       include "mpif.h"
 922       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 923 !      include "COMMON.MPI"
 924 #endif
 925 !      include "COMMON.CONTROL"
 926 !      include "COMMON.CHAIN"
 927 !      include "COMMON.IOUNITS"
 928 !      include "COMMON.PROTFILES"
 929 !      include "COMMON.NAMES"
 930 !      include "COMMON.VAR"
 931 !      include "COMMON.SBRIDGE"
 932 !      include "COMMON.GEO"
 933 !      include "COMMON.FFIELD"
 934 !      include "COMMON.ENEPS"
 935 !      include "COMMON.LOCAL"
 936 !      include "COMMON.WEIGHTS"
 937 !      include "COMMON.INTERACT"
 938 !      include "COMMON.FREE"
 939 !      include "COMMON.ENERGIES"
 940 !      include "COMMON.COMPAR"
 941 !      include "COMMON.PROT"
 942 !      include "COMMON.CONTACTS1"
 943       character(len=64) :: nazwa
 944       character(len=80) :: bxname,cxname
 945       character(len=64) :: bprotfile_temp
 946       character(len=3) :: liczba,licz
 947       character(len=2) :: licz2
 948       integer :: i,itj,ii,iii,j,k,l
 949       integer :: ixdrf,iret
 950       integer :: iscor,islice
 951       real(kind=8) :: rmsdev,efree,eini
 952       real(kind=4) :: csingle(3,nres*2)
 953       real(kind=8) :: energ
 954 !      integer ilen,iroof
 955 !      external ilen,iroof
 956       integer :: ir,ib,iparm
 957       integer :: isecstr(nres)
 958       write (licz2,'(bz,i2.2)') islice
 959       call opentmp(islice,ientout,bprotfile_temp)
 960       write (iout,*) "bprotfile_temp ",bprotfile_temp
 961       call flush(iout)
 962       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0 &
 963          .and. ensembles.eq.0) then
 964         close(ientout,status="delete")
 965         return
 966       endif
 967 #ifdef MPI
 968       write (liczba,'(bz,i3.3)') me
 969       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 970         if (.not.separate_parset) then
 971           bxname = prefix(:ilen(prefix))//liczba//".bx"
 972         else
 973           write (licz,'(bz,i3.3)') myparm
 974           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 975         endif
 976         open (ientin,file=bxname,status="unknown",&
 977           form="unformatted",access="direct",recl=lenrec1)
 978       endif
 979 #else
 980       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 981         if (nslice.eq.1) then
 982           bxname = prefix(:ilen(prefix))//".bx"
 983         else
 984           bxname = prefix(:ilen(prefix))// &
 985                  "_slice_"//licz2//".bx"
 986         endif
 987         open (ientin,file=bxname,status="unknown",&
 988           form="unformatted",access="direct",recl=lenrec1)
 989         write (iout,*) "Calculating energies; writing geometry",&
 990        " and energy components to ",bxname(:ilen(bxname))
 991       endif
 992 #if (defined(AIX) && !defined(JUBL))
 993         call xdrfopen_(ixdrf,cxname, "w", iret)
 994 #else
 995         call xdrfopen(ixdrf,cxname, "w", iret)
 996 #endif
 997         if (iret.eq.0) then
 998           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 999           cxfile=.false.
1000         endif
1001 !el      endif
1002 #endif
1003       if (indpdb.gt.0) then
1004         if (nslice.eq.1) then
1005 #ifdef MPI
1006          if (.not.separate_parset) then
1007            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot)) &
1008              //liczba//'.stat'
1009          else
1010            write (licz,'(bz,i3.3)') myparm
1011            statname=prefix(:ilen(prefix))//'_par'//licz//'_'// &
1012             pot(:ilen(pot))//liczba//'.stat'
1013          endif
1014
1015 #else
1016           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
1017 #endif
1018         else
1019 #ifdef MPI
1020          if (.not.separate_parset) then
1021           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))// &
1022             "_slice_"//licz2//liczba//'.stat'
1023          else
1024           write (licz,'(bz,i3.3)') myparm
1025           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))// &
1026             '_par'//licz//"_slice_"//licz2//liczba//'.stat'
1027          endif
1028 #else
1029           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot)) &
1030             //"_slice_"//licz2//'.stat'
1031 #endif
1032         endif
1033         open(istat,file=statname,status="unknown")
1034       endif
1035
1036 #ifdef MPI
1037       do i=1,scount(me)
1038 #else
1039       do i=1,ntot(islice)
1040 #endif
1041         read(ientout,rec=i,err=101) &
1042           ((csingle(l,k),l=1,3),k=1,nres),&
1043           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1044           nss,(ihpb(k),jhpb(k),k=1,nss),&
1045           eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
1046 !        write (iout,*) iR,ib,iparm,eini,efree
1047         do j=1,2*nres
1048           do k=1,3
1049             c(k,j)=csingle(k,j)
1050           enddo
1051         enddo
1052         call int_from_cart1(.false.)
1053         iscore=0
1054 !        write (iout,*) "Calling conf_compar",i
1055 !        call flush(iout)
1056          anatemp= 1.0d0/(beta_h(ib,iparm)*1.987D-3)
1057         if (indpdb.gt.0) then
1058           call conf_compar(i,.false.,.true.)
1059 !        else
1060 !            call elecont(.false.,ncont,icont,nnt,nct)
1061 !            call secondary2(.false.,.false.,ncont,icont,isecstr)
1062         endif
1063 !        write (iout,*) "Exit conf_compar",i
1064 !        call flush(iout)
1065         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i) &
1066           ((csingle(l,k),l=1,3),k=1,nres),&
1067           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1068           nss,(ihpb(k),jhpb(k),k=1,nss),&
1069 !     &    potE(i,iparm),-entfac(i),rms_nat,iscore
1070           potE(i,nparmset),-entfac(i),rms_nat,iscore
1071 !        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
1072 #ifndef MPI
1073         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),&
1074           -entfac(i),rms_nat,iscore)
1075 #endif
1076       enddo
1077       close(ientout,status="delete")
1078       close(istat)
1079       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
1080 #ifdef MPI
1081       call MPI_Barrier(WHAM_COMM,IERROR)
1082       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile &
1083          .and. ensembles.eq.0) return
1084       write (iout,*)
1085       if (bxfile .or. ensembles.gt.0) then
1086         if (nslice.eq.1) then
1087           if (.not.separate_parset) then
1088             bxname = prefix(:ilen(prefix))//".bx"
1089           else
1090             write (licz,'(bz,i3.3)') myparm
1091             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
1092           endif
1093         else
1094           if (.not.separate_parset) then
1095             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1096           else
1097             write (licz,'(bz,i3.3)') myparm
1098             bxname = prefix(:ilen(prefix))//"par_"//licz// &
1099               "_slice_"//licz2//".bx"
1100           endif
1101         endif
1102         open (ientout,file=bxname,status="unknown",&
1103             form="unformatted",access="direct",recl=lenrec1)
1104         write (iout,*) "Master is creating binary database ",&
1105          bxname(:ilen(bxname))
1106       endif
1107       if (cxfile) then
1108         if (nslice.eq.1) then
1109           if (.not.separate_parset) then
1110             cxname = prefix(:ilen(prefix))//".cx"
1111           else
1112             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
1113           endif
1114         else
1115           if (.not.separate_parset) then
1116             cxname = prefix(:ilen(prefix))// &
1117                    "_slice_"//licz2//".cx"
1118           else
1119             cxname = prefix(:ilen(prefix))//"_par"//licz// &
1120                    "_slice_"//licz2//".cx"
1121           endif
1122         endif
1123 #if (defined(AIX) && !defined(JUBL))
1124         call xdrfopen_(ixdrf,cxname, "w", iret)
1125 #else
1126         call xdrfopen(ixdrf,cxname, "w", iret)
1127 #endif
1128         if (iret.eq.0) then
1129           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
1130           cxfile=.false.
1131         endif
1132       endif
1133       do j=0,nprocs-1
1134         write (liczba,'(bz,i3.3)') j
1135         if (separate_parset) then
1136           write (licz,'(bz,i3.3)') myparm
1137           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
1138         else
1139           bxname = prefix(:ilen(prefix))//liczba//".bx"
1140         endif
1141         open (ientin,file=bxname,status="unknown",&
1142           form="unformatted",access="direct",recl=lenrec1)
1143         write (iout,*) "Master is reading conformations from ",&
1144          bxname(:ilen(bxname))
1145         iii = 0
1146 !        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
1147 !        call flush(iout)
1148         do i=indstart(j),indend(j)
1149           iii = iii+1
1150           read(ientin,rec=iii,err=101) &
1151             ((csingle(l,k),l=1,3),k=1,nres),&
1152             ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1153             nss,(ihpb(k),jhpb(k),k=1,nss),&
1154             eini,efree,rmsdev,iscor
1155           if (bxfile .or. ensembles.gt.0) then
1156             write (ientout,rec=i) &
1157               ((csingle(l,k),l=1,3),k=1,nres),&
1158               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1159               nss,(ihpb(k),jhpb(k),k=1,nss),&
1160               eini,efree,rmsdev,iscor
1161           endif
1162           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
1163 #ifdef DEBUG
1164           do k=1,2*nres
1165             do l=1,3
1166               c(l,k)=csingle(l,k)
1167             enddo
1168           enddo
1169           call int_from_cart1(.false.)
1170           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
1171           write (iout,*) "The Cartesian geometry is:"
1172           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1173           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1174           write (iout,*) "The internal geometry is:"
1175           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1176           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1177           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1178           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1179           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1180           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1181           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
1182           write (iout,'(f10.5,i5)') rmsdev,iscor
1183 #endif
1184         enddo ! i
1185         write (iout,*) iii," conformations (from",indstart(j)," to",&
1186          indend(j),") read from ",&
1187          bxname(:ilen(bxname))
1188         close (ientin,status="delete")
1189       enddo ! j
1190       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
1191 #if (defined(AIX) && !defined(JUBL))
1192       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
1193 #else
1194       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
1195 #endif
1196 #endif
1197       return
1198   101 write (iout,*) "Error in scratchfile."
1199       call flush(iout)
1200       return 1
1201       end subroutine write_dbase
1202 !-------------------------------------------------------------------------------
1203       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
1204 !      implicit none
1205 !      include "DIMENSIONS"
1206 !      include "DIMENSIONS.ZSCOPT"
1207 !      include "DIMENSIONS.FREE"
1208 !      include "DIMENSIONS.COMPAR"
1209 #ifdef MPI
1210       include "mpif.h"
1211       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1212 !      include "COMMON.MPI"
1213 #endif
1214 !      include "COMMON.CONTROL"
1215 !      include "COMMON.CHAIN"
1216 !      include "COMMON.IOUNITS"
1217 !      include "COMMON.PROTFILES"
1218 !      include "COMMON.NAMES"
1219 !      include "COMMON.VAR"
1220 !      include "COMMON.SBRIDGE"
1221 !      include "COMMON.GEO"
1222 !      include "COMMON.FFIELD"
1223 !      include "COMMON.ENEPS"
1224 !      include "COMMON.LOCAL"
1225 !      include "COMMON.WEIGHTS"
1226 !      include "COMMON.INTERACT"
1227 !      include "COMMON.FREE"
1228 !      include "COMMON.ENERGIES"
1229 !      include "COMMON.COMPAR"
1230 !      include "COMMON.PROT"
1231       integer :: i,j,itmp,iscor,iret,ixdrf
1232       real(kind=8) :: rmsdev,efree,eini
1233       real(kind=4) :: csingle(3,nres*2),xoord(3,2*nres+2)
1234       real(kind=4) :: prec
1235
1236 !      write (iout,*) "cxwrite"
1237 !      call flush(iout)
1238       prec=10000.0
1239       do i=1,nres
1240        do j=1,3
1241         xoord(j,i)=csingle(j,i)
1242        enddo
1243       enddo
1244       do i=nnt,nct
1245        do j=1,3
1246         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
1247        enddo
1248       enddo
1249
1250       itmp=nres+nct-nnt+1
1251
1252 !      write (iout,*) "itmp",itmp
1253 !      call flush(iout)
1254 #if (defined(AIX) && !defined(JUBL))
1255       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
1256
1257 !      write (iout,*) "xdrf3dfcoord"
1258 !      call flush(iout)
1259       call xdrfint_(ixdrf, nss, iret)
1260       do j=1,nss
1261         call xdrfint_(ixdrf, ihpb(j), iret)
1262         call xdrfint_(ixdrf, jhpb(j), iret)
1263       enddo
1264       call xdrffloat_(ixdrf,real(eini),iret)
1265       call xdrffloat_(ixdrf,real(efree),iret)
1266       call xdrffloat_(ixdrf,real(rmsdev),iret)
1267       call xdrfint_(ixdrf,iscor,iret)
1268 #else
1269       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
1270
1271       call xdrfint(ixdrf, nss, iret)
1272       do j=1,nss
1273         call xdrfint(ixdrf, ihpb(j), iret)
1274         call xdrfint(ixdrf, jhpb(j), iret)
1275       enddo
1276       call xdrffloat(ixdrf,real(eini),iret)
1277       call xdrffloat(ixdrf,real(efree),iret)
1278       call xdrffloat(ixdrf,real(rmsdev),iret)
1279       call xdrfint(ixdrf,iscor,iret)
1280 #endif
1281
1282       return
1283       end subroutine cxwrite
1284 !-------------------------------------------------------------------------------
1285 ! slices.F
1286 !-------------------------------------------------------------------------------
1287       subroutine set_slices(is,ie,ts,te,iR,ib,iparm)
1288 !      implicit none
1289 !      include 'DIMENSIONS'
1290 !      include 'DIMENSIONS.ZSCOPT'
1291 !      include 'DIMENSIONS.FREE'
1292 !      include 'COMMON.IOUNITS'
1293 !      include 'COMMON.PROTFILES'
1294 !      include 'COMMON.OBCINKA'
1295 !      include 'COMMON.PROT'
1296       integer :: islice,iR,ib,iparm
1297       integer :: is(MaxSlice),ie(MaxSlice),nrec_slice
1298       real(kind=8) :: ts(MaxSlice),te(MaxSlice),time_slice
1299
1300       do islice=1,nslice
1301         if (time_end_collect(iR,ib,iparm).ge.1.0d10) then
1302           ts(islice)=time_start_collect(iR,ib,iparm)
1303           te(islice)=time_end_collect(iR,ib,iparm)
1304           nrec_slice=(rec_end(iR,ib,iparm)- &
1305              rec_start(iR,ib,iparm)+1)/nslice
1306           is(islice)=rec_start(iR,ib,iparm)+(islice-1)*nrec_slice
1307           ie(islice)=rec_start(iR,ib,iparm)+islice*nrec_slice-1
1308         else
1309           time_slice=(time_end_collect(iR,ib,iparm) &
1310           -time_start_collect(iR,ib,iparm))/nslice
1311           ts(islice)=time_start_collect(iR,ib,iparm)+(islice-1)* &
1312            time_slice
1313           te(islice)=time_start_collect(iR,ib,iparm)+islice*time_slice
1314           is(islice)=rec_start(iR,ib,iparm)
1315           ie(islice)=rec_end(iR,ib,iparm)
1316         endif
1317       enddo
1318
1319       write (iout,*) "nrec_slice",nrec_slice," time_slice",time_slice
1320       write (iout,*) "is",(is(islice),islice=1,nslice)
1321       write (iout,*) "ie",(ie(islice),islice=1,nslice)
1322       write (iout,*) "rec_start",&
1323         rec_start(iR,ib,iparm)," rec_end",rec_end(iR,ib,iparm)
1324       write (iout,*) "ts",(ts(islice),islice=1,nslice)
1325       write (iout,*) "te",(te(islice),islice=1,nslice)
1326       write (iout,*) "time_start",&
1327         time_start_collect(iR,ib,iparm)," time_end",&
1328         time_end_collect(iR,ib,iparm)
1329       call flush(iout)
1330
1331       return
1332       end subroutine set_slices
1333 !-----------------------------------------------------------------------------
1334       integer function slice(irecord,time,is,ie,ts,te)
1335 !      implicit none
1336 !      include 'DIMENSIONS'
1337 !      include 'DIMENSIONS.ZSCOPT'
1338 !      include 'DIMENSIONS.FREE'
1339 !      include 'COMMON.IOUNITS'
1340 !      include 'COMMON.PROTFILES'
1341 !      include 'COMMON.OBCINKA'
1342 !      include 'COMMON.PROT'
1343       integer :: is(MaxSlice),ie(MaxSlice),nrec_slice
1344       real(kind=8) :: ts(MaxSlice),te(MaxSlice),time_slice
1345       integer :: i,ii,irecord
1346       real(kind=8) :: time
1347
1348 !      write (iout,*) "within slice nslice",nslice
1349 !      call flush(iout)
1350       if (irecord.lt.is(1) .or. time.lt.ts(1)) then
1351         ii=0
1352       else
1353         ii=1
1354         do while (ii.le.nslice .and. &
1355                  (irecord.lt.is(ii) .or. irecord.gt.ie(ii) .or. &
1356                  time.lt.ts(ii) .or. time.gt.te(ii)) )
1357 !          write (iout,*) "ii",ii,time,ts(ii)
1358 !          call flush(iout)
1359           ii=ii+1
1360         enddo
1361       endif
1362 !      write (iout,*) "end: ii",ii
1363 !      call flush(iout)
1364       slice=ii
1365       return
1366       end function slice
1367 !-----------------------------------------------------------------------------
1368 ! enecalc1.F
1369 !-----------------------------------------------------------------------------
1370       logical function conf_check(ii,iprint)
1371
1372       use names, only:ntyp1
1373       use geometry_data
1374       use energy_data, only:itype,dsc
1375       use geometry, only:int_from_cart1
1376 !      use
1377 !      include "DIMENSIONS"
1378 !      include "DIMENSIONS.ZSCOPT"
1379 !      include "DIMENSIONS.FREE"
1380 !#ifdef MPI
1381 !      use MPI_data
1382 !      include "mpif.h"
1383 !      include "COMMON.MPI"
1384 !#endif
1385 !      include "COMMON.CHAIN"
1386 !      include "COMMON.IOUNITS"
1387 !      include "COMMON.PROTFILES"
1388 !      include "COMMON.NAMES"
1389 !      include "COMMON.VAR"
1390 !      include "COMMON.SBRIDGE"
1391 !      include "COMMON.GEO"
1392 !      include "COMMON.FFIELD"
1393 !      include "COMMON.ENEPS"
1394 !      include "COMMON.LOCAL"
1395 !      include "COMMON.WEIGHTS"
1396 !      include "COMMON.INTERACT"
1397 !      include "COMMON.FREE"
1398 !      include "COMMON.ENERGIES"
1399 !      include "COMMON.CONTROL"
1400 !      include "COMMON.TORCNSTR"
1401 !      implicit none
1402 #ifdef MPI
1403       include "mpif.h"
1404       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1405 #endif
1406       integer :: j,k,l,ii,itj,iprint
1407       if (.not. check_conf) then
1408         conf_check=.true.
1409         return
1410       endif
1411       call int_from_cart1(.false.)
1412       do j=nnt+1,nct
1413         if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. &
1414           (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
1415           if (iprint.gt.0) &
1416           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
1417             " for conformation",ii
1418           if (iprint.gt.1) then
1419             write (iout,*) "The Cartesian geometry is:"
1420             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1421             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1422             write (iout,*) "The internal geometry is:"
1423             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1424             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1425             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1426             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1427             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1428             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1429           endif
1430           if (iprint.gt.0) write (iout,*) &
1431             "This conformation WILL NOT be added to the database."
1432           conf_check=.false.
1433           return
1434         endif
1435       enddo
1436       do j=nnt,nct
1437         itj=itype(j)
1438         if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. &
1439            (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
1440           if (iprint.gt.0) &
1441           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
1442            " for conformation",ii
1443           if (iprint.gt.1) then
1444             write (iout,*) "The Cartesian geometry is:"
1445             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1446             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1447             write (iout,*) "The internal geometry is:"
1448             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1449             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1450             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1451             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1452             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1453             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1454           endif
1455           if (iprint.gt.0) write (iout,*) &
1456             "This conformation WILL NOT be added to the database."
1457           conf_check=.false.
1458           return
1459         endif
1460       enddo
1461       do j=3,nres
1462         if (theta(j).le.0.0d0) then
1463           if (iprint.gt.0) &
1464           write (iout,*) "Zero theta angle(s) in conformation",ii
1465           if (iprint.gt.1) then
1466             write (iout,*) "The Cartesian geometry is:"
1467             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1468             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1469             write (iout,*) "The internal geometry is:"
1470             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1471             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1472             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1473             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1474             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1475             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1476           endif
1477           if (iprint.gt.0) write (iout,*) &
1478             "This conformation WILL NOT be added to the database."
1479           conf_check=.false.
1480           return
1481         endif
1482         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
1483       enddo
1484       conf_check=.true.
1485 !      write (iout,*) "conf_check passed",ii
1486       return
1487       end function conf_check
1488 !-----------------------------------------------------------------------------
1489       end module io_database