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