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       print *,"after xdrf3dcoord", iret
 547 #endif
 548 #ifdef DEBUG
 549       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,2*nres+2)
 550 #endif
 551 !#undef DEBUG
 552       if (iret.eq.0) exit
 553       if (itmp .ne. nres + nct - nnt + 1) then
 554         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
 555         call flush(iout)
 556         exit
 557       endif
 558
 559       time=rtime
 560 !      write (iout,*) "calling slice" !d
 561       call flush(iout) !d
 562       islice=slice(nstep(itraj),time,is,ie,ts,te)
 563 !      write (iout,*) "islice",islice !d
 564       call flush(iout) !d
 565
 566       do i=1,nres
 567         do j=1,3
 568           c(j,i)=xoord(j,i)
 569         enddo
 570       enddo
 571       do i=1,nct-nnt+1
 572         do j=1,3
 573           c(j,i+nres+nnt-1)=xoord(j,i+nres)
 574         enddo
 575       enddo
 576
 577       if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset &
 578           .or. iset.eq.myparm)) then
 579         ii=ii+1
 580         kk(islice)=kk(islice)+1
 581         mm(islice)=mm(islice)+1
 582         if (mod(nstep(itraj),isampl(iparm)).eq.0 .and. &
 583            conf_check(ll(islice)+1,1)) then
 584           if (replica(iparm)) then
 585              rt_bath=1.0d0/(rt_bath*1.987D-3)
 586              do i=1,nT_h(iparm)
 587                if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
 588                  iib = i
 589                  goto 22
 590                endif
 591              enddo
 592   22         continue
 593              if (i.gt.nT_h(iparm)) then
 594                write (iout,*) "Error - temperature of conformation",&
 595                ii,1.0d0/(rt_bath*1.987D-3),&
 596                " does not match any of the list"
 597                write (iout,*) &
 598                 1.0d0/(rt_bath*1.987D-3),&
 599                 (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 600                call flush(iout)
 601 !               exit
 602 !               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 603                ii=ii-1
 604                kk(islice)=kk(islice)-1
 605                mm(islice)=mm(islice)-1
 606                goto 112
 607              endif
 608           else
 609             iib = ib
 610           endif
 611
 612           efree=0.0d0
 613           jj(islice)=jj(islice)+1
 614           if (umbrella(iparm)) then
 615             snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 616           else if (hamil_rep) then
 617             snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
 618           else
 619             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 620           endif
 621           ll(islice)=ll(islice)+1
 622 #ifdef DEBUG
 623           write (iout,*) "Writing conformation, record",ll(islice)
 624           write (iout,*) "ib",ib," iib",iib
 625           write (iout,*) "ntraj",ntraj," itraj",itraj,&
 626             " nstep",nstep(itraj)
 627           write (iout,*) "pote",rpotE," time",rtime
 628 !          if (replica(iparm)) then
 629 !            write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
 630 !            write (iout,*) "TEMP list"
 631 !            write (iout,*)
 632 !     &       (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 633 !          endif
 634           write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 635 !          write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 636 !          write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 637           call flush(iout)
 638 #endif
 639           if (islice.ne.islice1) then
 640 !            write (iout,*) "islice",islice," islice1",islice1
 641             close(ientout)
 642 !            write (iout,*) "Closing file ",
 643 !     &          bprotfile_temp(:ilen(bprotfile_temp))
 644             call opentmp(islice,ientout,bprotfile_temp)
 645 !            write (iout,*) "Opening file ",
 646 !     &          bprotfile_temp(:ilen(bprotfile_temp))
 647             islice1=islice
 648           endif
 649           if (umbrella(iparm)) then
 650             write(ientout,rec=ll(islice)) &
 651               ((xoord(l,k),l=1,3),k=1,nres),&
 652               ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),&
 653               nss,(ihpb(k),jhpb(k),k=1,nss),&
 654               rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),&
 655               iset,iib,iparm
 656           else if (hamil_rep) then
 657             write(ientout,rec=ll(islice)) &
 658               ((xoord(l,k),l=1,3),k=1,nres),&
 659               ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),&
 660               nss,(ihpb(k),jhpb(k),k=1,nss),&
 661               rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),&
 662               iR,iib,iset
 663           else
 664             write(ientout,rec=ll(islice)) &
 665               ((xoord(l,k),l=1,3),k=1,nres),&
 666               ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),&
 667               nss,(ihpb(k),jhpb(k),k=1,nss),&
 668               rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),&
 669               iR,iib,iparm
 670           endif
 671 #ifdef DEBUG
 672           call int_from_cart1(.false.)
 673           write (iout,*) "Writing conformation, record",ll(islice)
 674           write (iout,*) "Cartesian coordinates"
 675           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 676           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 677           write (iout,*) "Internal coordinates"
 678           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 679           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 680           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 681           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 682           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 683           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 684           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 685 !          write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
 686           write (iout,'(16i5)') iscor
 687           call flush(iout)
 688 #endif
 689         endif
 690       endif
 691
 692   112 continue
 693
 694       enddo
 695       close(ientout)
 696 #if (defined(AIX) && !defined(JUBL))
 697       call xdrfclose_(ixdrf, iret)
 698 #else
 699       call xdrfclose(ixdrf, iret)
 700 #endif
 701       write (iout,'(i10," trajectories found in file.")') ntraj+1
 702       write (iout,'(a)') "Numbers of steps in trajectories:"
 703       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 704       write (iout,*) ii," conformations read from file",&
 705          nazwa(:ilen(nazwa))
 706       do islice=1,nslice
 707         write (iout,*) mm(islice)," conformations read so far, slice",&
 708           islice
 709         write (iout,*) ll(islice),&
 710         " conformations stored so far, slice",islice
 711       enddo
 712       call flush(iout)
 713       print *,"before cxread return"
 714 !#undef DEBUG
 715       return
 716       end subroutine cxread
 717 !--------------------------------------------------------------------------------
 718 ! xread.F
 719 !--------------------------------------------------------------------------------
 720       subroutine xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
 721
 722       use geometry_data
 723 !      implicit none
 724 !      include "DIMENSIONS"
 725 !      include "DIMENSIONS.ZSCOPT"
 726 !      include "DIMENSIONS.FREE"
 727       use MPI_data, only:nprocs
 728 #ifdef MPI
 729       include "mpif.h"
 730       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 731 !      include "COMMON.MPI"
 732 #endif
 733       integer,parameter :: MaxTraj=2050
 734 !      include "COMMON.CHAIN"
 735 !      include "COMMON.IOUNITS"
 736 !      include "COMMON.PROTFILES"
 737 !      include "COMMON.NAMES"
 738 !      include "COMMON.VAR"
 739 !      include "COMMON.GEO"
 740 !      include "COMMON.ENEPS"
 741 !      include "COMMON.PROT"
 742 !      include "COMMON.INTERACT"
 743 !      include "COMMON.FREE"
 744 !      include "COMMON.SBRIDGE"
 745 !      include "COMMON.OBCINKA"
 746       real(kind=4) :: csingle(3,nres*2)
 747       character(len=64) :: nazwa,bprotfile_temp
 748       integer :: i,j,k,l,ii,jj(nslice),kk(nslice),ll(nslice),&
 749         mm(nslice) !(maxslice)
 750       integer :: iscor,islice,islice1 !el,slice
 751       real(kind=8) :: energ
 752 !      integer ilen,iroof
 753 !      external ilen,iroof
 754       real(kind=8) :: rmsdev,energia(0:n_ene),efree,eini,temp
 755 !el      real(kind=8) :: rmsdev,energia(0:max_eneW),efree,eini,temp
 756       real(kind=8) :: prop(nQ) !(maxQ)
 757       integer :: ntot_all(0:nprocs-1)!(0:maxprocs-1)
 758       integer :: iparm,ib,iib,ir,nprop,nthr
 759       real(kind=8) :: etot,time,ts(nslice),te(nslice)
 760       integer :: is(nslice),ie(nslice),itraj,ntraj,it,iset
 761       integer :: nstep(0:MaxTraj-1)
 762       logical :: lerr
 763
 764       call set_slices(is,ie,ts,te,iR,ib,iparm)
 765       do i=1,nQ
 766         prop(i)=0.0d0
 767       enddo
 768       do i=0,MaxTraj-1
 769         nstep(i)=0
 770       enddo
 771       ntraj=0
 772       it=0
 773       islice1=1
 774       call opentmp(islice1,ientout,bprotfile_temp)
 775       do while (.true.)
 776         if (replica(iparm)) then
 777           if (hamil_rep .or. umbrella(iparm)) then
 778           read (ientin,*,end=1112,err=1112) time,eini,&
 779             etot,temp,nss,(ihpb(j),jhpb(j),j=1,nss),&
 780             nprop,(prop(j),j=1,nprop),iset
 781           else
 782           read (ientin,*,end=1112,err=1112) time,eini,&
 783             etot,temp,nss,(ihpb(j),jhpb(j),j=1,nss),&
 784             nprop,(prop(j),j=1,nprop)
 785           endif
 786           temp=1.0d0/(temp*1.987D-3)
 787 !           write (iout,*) time,eini,etot,nss,
 788 !     &     (ihpb(j),jhpb(j),j=1,nss),(prop(j),j=1,nprop)
 789 !           call flush(iout)
 790            do i=1,nT_h(iparm)
 791              if (beta_h(i,iparm).eq.temp) then
 792                iib = i
 793                goto 22
 794              endif
 795            enddo
 796   22       continue
 797            if (i.gt.nT_h(iparm)) then
 798              write (iout,*) "Error - temperature of conformation",&
 799              ii,1.0d0/(temp*1.987D-3),&
 800              " does not match any of the list"
 801              write (iout,*) &
 802               1.0d0/(temp*1.987D-3),&
 803               (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 804              call flush(iout)
 805 #ifdef MPI
 806              call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 807 #endif
 808            endif
 809         else
 810            read (ientin,*,end=1112,err=1112) time,eini,&
 811              etot,nss,(ihpb(j),jhpb(j),j=1,nss),&
 812              nprop,(prop(j),j=1,nprop)
 813              iib = ib
 814         endif
 815         itraj=mod(it,totraj(iR,iparm))
 816 !        write (*,*) "ii",ii," itraj",itraj
 817 !        call flush(iout)
 818         it=it+1
 819         if (itraj.gt.ntraj) ntraj=itraj
 820         nstep(itraj)=nstep(itraj)+1
 821         islice=slice(nstep(itraj),time,is,ie,ts,te)
 822         read (ientin,'(8f10.5)',end=1112,err=1112) &
 823           ((csingle(l,k),l=1,3),k=1,nres),&
 824           ((csingle(l,k+nres),l=1,3),k=nnt,nct)
 825         efree=0.0d0
 826         if (islice.gt.0 .and. islice.le.nslice) then
 827         ii=ii+1
 828         kk(islice)=kk(islice)+1
 829         mm(islice)=mm(islice)+1
 830         if (mod(nstep(itraj),isampl(iparm)).eq.0) then
 831         jj(islice)=jj(islice)+1
 832         if (hamil_rep) then
 833           snk(iR,iib,iset,islice)=snk(iR,iib,iset,islice)+1
 834         else if (umbrella(iparm)) then
 835           snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 836         else
 837           snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 838         endif
 839         ll(islice)=ll(islice)+1
 840 !         write (iout,*) ii,kk,jj,ll,eini,(prop(j),j=1,nprop)
 841 #ifdef DEBUG
 842 !        write (iout,*) "Writing conformation, record",ll(islice)
 843 !        write (iout,*) "ib",ib," iib",iib
 844          if (replica(iparm)) then
 845            write (iout,*) "TEMP",1.0d0/(temp*1.987D-3)
 846            write (iout,*) "TEMP list"
 847            write (iout,*) &
 848             (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 849          endif
 850          call flush(iout)
 851 #endif
 852 !         write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 853 !         write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 854 !         write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 855 !         call flush(iout)
 856          if (islice.ne.islice1) then
 857 !             write (iout,*) "islice",islice," islice1",islice1
 858              close(ientout)
 859 !             write (iout,*) "Closing file ",
 860 !     &             bprotfile_temp(:ilen(bprotfile_temp))
 861              call opentmp(islice,ientout,bprotfile_temp)
 862 !             write (iout,*) "Opening file ",
 863 !     &             bprotfile_temp(:ilen(bprotfile_temp))
 864 !             call flush(iout)
 865              islice1=islice
 866          endif
 867          write(ientout,rec=ll(islice)) &
 868            ((csingle(l,k),l=1,3),k=1,nres),&
 869            ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
 870            nss,(ihpb(k),jhpb(k),k=1,nss),&
 871            eini,efree,rmsdev,(prop(i),i=1,nQ),iR,iib,iparm
 872 #ifdef DEBUG
 873          do i=1,2*nres
 874            do j=1,3
 875              c(j,i)=csingle(j,i)
 876            enddo
 877          enddo
 878          call int_from_cart1(.false.)
 879          write (iout,*) "Writing conformation, record",ll(islice)
 880          write (iout,*) "Cartesian coordinates"
 881          write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 882          write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 883          write (iout,*) "Internal coordinates"
 884          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 885          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 886          write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 887          write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 888          write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 889          write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 890          write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 891 !         write (iout,'(8f10.5)') (prop(j),j=1,nQ)
 892          write (iout,'(16i5)') iscor
 893          call flush(iout)
 894 #endif
 895          endif
 896          endif
 897        enddo
 898  1112  continue
 899        close(ientout)
 900        write (iout,'(i10," trajectories found in file.")') ntraj+1
 901        write (iout,'(a)') "Numbers of steps in trajectories:"
 902        write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 903        write (iout,*) ii," conformations read from file",&
 904          nazwa(:ilen(nazwa))
 905        write (iout,*) mm(islice)," conformations read so far, slice",&
 906           islice
 907        write (iout,*) ll(islice)," conformations stored so far, slice",&
 908          islice
 909        call flush(iout)
 910        return
 911        end subroutine xread
 912 !--------------------------------------------------------------------------------
 913 ! enecalc1.F
 914 !--------------------------------------------------------------------------------
 915       subroutine write_dbase(islice,*)
 916
 917       use geometry_data
 918       use control_data, only:indpdb
 919       use w_compar_data
 920       use conform_compar, only:conf_compar,rmsnat,qwolynes
 921       use energy_data, only:icont,ncont,nnt,nct,maxcont!,&
 922 !      implicit none
 923 !      include "DIMENSIONS"
 924 !      include "DIMENSIONS.ZSCOPT"
 925 !      include "DIMENSIONS.FREE"
 926 !      include "DIMENSIONS.COMPAR"
 927       use geometry, only:int_from_cart1
 928 #ifdef MPI
 929       include "mpif.h"
 930       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 931 !      include "COMMON.MPI"
 932 #endif
 933 !      include "COMMON.CONTROL"
 934 !      include "COMMON.CHAIN"
 935 !      include "COMMON.IOUNITS"
 936 !      include "COMMON.PROTFILES"
 937 !      include "COMMON.NAMES"
 938 !      include "COMMON.VAR"
 939 !      include "COMMON.SBRIDGE"
 940 !      include "COMMON.GEO"
 941 !      include "COMMON.FFIELD"
 942 !      include "COMMON.ENEPS"
 943 !      include "COMMON.LOCAL"
 944 !      include "COMMON.WEIGHTS"
 945 !      include "COMMON.INTERACT"
 946 !      include "COMMON.FREE"
 947 !      include "COMMON.ENERGIES"
 948 !      include "COMMON.COMPAR"
 949 !      include "COMMON.PROT"
 950 !      include "COMMON.CONTACTS1"
 951       character(len=64) :: nazwa
 952       character(len=80) :: bxname,cxname
 953       character(len=64) :: bprotfile_temp
 954       character(len=3) :: liczba,licz
 955       character(len=2) :: licz2
 956       integer :: i,itj,ii,iii,j,k,l
 957       integer :: ixdrf,iret
 958       integer :: iscor,islice
 959       real(kind=8) :: rmsdev,efree,eini,qnat2
 960       real(kind=4) :: csingle(3,nres*2)
 961       real(kind=8) :: energ
 962
 963 !      integer ilen,iroof
 964 !      external ilen,iroof
 965       integer :: ir,ib,iparm
 966       integer :: isecstr(nres)
 967       logical :: test
 968       write (licz2,'(bz,i2.2)') islice
 969       call opentmp(islice,ientout,bprotfile_temp)
 970       write (iout,*) "bprotfile_temp ",bprotfile_temp
 971       call flush(iout)
 972       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0 &
 973          .and. ensembles.eq.0) then
 974         close(ientout,status="delete")
 975         return
 976       endif
 977 #ifdef MPI
 978       write (liczba,'(bz,i3.3)') me
 979       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 980         if (.not.separate_parset) then
 981           bxname = prefix(:ilen(prefix))//liczba//".bx"
 982         else
 983           write (licz,'(bz,i3.3)') myparm
 984           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 985         endif
 986         print *,bxname
 987         open (ientin,file=bxname,status="unknown",&
 988           form="unformatted",access="direct",recl=lenrec1)
 989       endif
 990 #else
 991       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 992         if (nslice.eq.1) then
 993           bxname = prefix(:ilen(prefix))//".bx"
 994         else
 995           bxname = prefix(:ilen(prefix))// &
 996                  "_slice_"//licz2//".bx"
 997         endif
 998         open (ientin,file=bxname,status="unknown",&
 999           form="unformatted",access="direct",recl=lenrec1)
1000         write (iout,*) "Calculating energies; writing geometry",&
1001        " and energy components to ",bxname(:ilen(bxname))
1002       endif
1003 #if (defined(AIX) && !defined(JUBL))
1004         call xdrfopen_(ixdrf,cxname, "w", iret)
1005 #else
1006         call xdrfopen(ixdrf,cxname, "w", iret)
1007 #endif
1008         if (iret.eq.0) then
1009           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
1010           cxfile=.false.
1011         endif
1012 !el      endif
1013 #endif
1014       print *,indpdb
1015       if (indpdb.gt.0) then
1016         if (nslice.eq.1) then
1017 #ifdef MPI
1018          if (.not.separate_parset) then
1019            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot)) &
1020              //liczba//'.stat'
1021          else
1022            write (licz,'(bz,i3.3)') myparm
1023            statname=prefix(:ilen(prefix))//'_par'//licz//'_'// &
1024             pot(:ilen(pot))//liczba//'.stat'
1025          endif
1026          print *,statname
1027 #else
1028           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
1029 #endif
1030         else
1031 #ifdef MPI
1032          if (.not.separate_parset) then
1033           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))// &
1034             "_slice_"//licz2//liczba//'.stat'
1035          else
1036           write (licz,'(bz,i3.3)') myparm
1037           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))// &
1038             '_par'//licz//"_slice_"//licz2//liczba//'.stat'
1039          endif
1040 #else
1041           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot)) &
1042             //"_slice_"//licz2//'.stat'
1043 #endif
1044         endif
1045         print *,istat,statname
1046         open(istat,file=statname,status="unknown")
1047       endif
1048       print *,"Tu dochodze"
1049       print *,scount(me)
1050 #ifdef MPI
1051       do i=1,scount(me)
1052 #else
1053       do i=1,ntot(islice)
1054 #endif
1055         print *,"before ientout read"
1056         read(ientout,rec=i,err=101) &
1057           ((csingle(l,k),l=1,3),k=1,nres),&
1058           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1059           nss,(ihpb(k),jhpb(k),k=1,nss),&
1060           eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
1061 !        write (iout,*) iR,ib,iparm,eini,efree
1062         do j=1,2*nres
1063           do k=1,3
1064             c(k,j)=csingle(k,j)
1065           enddo
1066         enddo
1067         call int_from_cart1(.false.)
1068         iscore=0
1069 !        write (iout,*) "Calling conf_compar",i
1070 !        call flush(iout)
1071          anatemp= 1.0d0/(beta_h(ib,iparm)*1.987D-3)
1072         print *,"before conf_compar"
1073         if (indpdb.gt.0) then
1074         print *,"just before conf_compar",i
1075         print *,icont,ncont,nnt,nct,"maxcont",maxcont
1076         test=.false.
1077 !          call conf_compar(i,.false.,.true.)
1078 !          call conf_compar(i)
1079 !           call rmsnat(i)
1080            rms_nat=rmsnat(i)
1081            qnat2=qwolynes(0,0)
1082          print *,"just after conf_compar"
1083 !        else
1084 !            call elecont(.false.,ncont,icont,nnt,nct)
1085 !            call secondary2(.false.,.false.,ncont,icont,isecstr)
1086         endif
1087 !        write (iout,*) "Exit conf_compar",i
1088 !        call flush(iout)
1089          print *,"before ientin"
1090         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i) &
1091           ((csingle(l,k),l=1,3),k=1,nres),&
1092           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1093           nss,(ihpb(k),jhpb(k),k=1,nss),&
1094 !     &    potE(i,iparm),-entfac(i),rms_nat,iscore
1095           potE(i,nparmset),-entfac(i),rms_nat,iscore
1096 !        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
1097 #ifndef MPI
1098         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),&
1099           -entfac(i),rms_nat,iscore)
1100 #endif
1101       enddo
1102       close(ientout,status="delete")
1103       close(istat)
1104       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
1105 #ifdef MPI
1106       print *,"before MPI_barrier"
1107       call MPI_Barrier(WHAM_COMM,IERROR)
1108       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile &
1109          .and. ensembles.eq.0) return
1110       write (iout,*)
1111       if (bxfile .or. ensembles.gt.0) then
1112         if (nslice.eq.1) then
1113           if (.not.separate_parset) then
1114             bxname = prefix(:ilen(prefix))//".bx"
1115           else
1116             write (licz,'(bz,i3.3)') myparm
1117             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
1118           endif
1119         else
1120           if (.not.separate_parset) then
1121             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1122           else
1123             write (licz,'(bz,i3.3)') myparm
1124             bxname = prefix(:ilen(prefix))//"par_"//licz// &
1125               "_slice_"//licz2//".bx"
1126           endif
1127         endif
1128         open (ientout,file=bxname,status="unknown",&
1129             form="unformatted",access="direct",recl=lenrec1)
1130         write (iout,*) "Master is creating binary database ",&
1131          bxname(:ilen(bxname))
1132       endif
1133       if (cxfile) then
1134         if (nslice.eq.1) then
1135           if (.not.separate_parset) then
1136             cxname = prefix(:ilen(prefix))//".cx"
1137           else
1138             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
1139           endif
1140         else
1141           if (.not.separate_parset) then
1142             cxname = prefix(:ilen(prefix))// &
1143                    "_slice_"//licz2//".cx"
1144           else
1145             cxname = prefix(:ilen(prefix))//"_par"//licz// &
1146                    "_slice_"//licz2//".cx"
1147           endif
1148         endif
1149 #if (defined(AIX) && !defined(JUBL))
1150         call xdrfopen_(ixdrf,cxname, "w", iret)
1151 #else
1152         call xdrfopen(ixdrf,cxname, "w", iret)
1153 #endif
1154         if (iret.eq.0) then
1155           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
1156           cxfile=.false.
1157         endif
1158       endif
1159       do j=0,nprocs-1
1160         write (liczba,'(bz,i3.3)') j
1161         if (separate_parset) then
1162           write (licz,'(bz,i3.3)') myparm
1163           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
1164         else
1165           bxname = prefix(:ilen(prefix))//liczba//".bx"
1166         endif
1167         open (ientin,file=bxname,status="unknown",&
1168           form="unformatted",access="direct",recl=lenrec1)
1169         write (iout,*) "Master is reading conformations from ",&
1170          bxname(:ilen(bxname))
1171         iii = 0
1172 !        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
1173 !        call flush(iout)
1174         do i=indstart(j),indend(j)
1175           iii = iii+1
1176           read(ientin,rec=iii,err=101) &
1177             ((csingle(l,k),l=1,3),k=1,nres),&
1178             ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1179             nss,(ihpb(k),jhpb(k),k=1,nss),&
1180             eini,efree,rmsdev,iscor
1181           if (bxfile .or. ensembles.gt.0) then
1182             write (ientout,rec=i) &
1183               ((csingle(l,k),l=1,3),k=1,nres),&
1184               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1185               nss,(ihpb(k),jhpb(k),k=1,nss),&
1186               eini,efree,rmsdev,iscor
1187           endif
1188 !           print *,"before cxwrite"
1189           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
1190 !          print *,"after cxwrite"
1191 #ifdef DEBUG
1192           do k=1,2*nres
1193             do l=1,3
1194               c(l,k)=csingle(l,k)
1195             enddo
1196           enddo
1197           call int_from_cart1(.false.)
1198           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
1199           write (iout,*) "The Cartesian geometry is:"
1200           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1201           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1202           write (iout,*) "The internal geometry is:"
1203           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1204           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1205           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1206           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1207           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1208           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1209           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
1210           write (iout,'(f10.5,i5)') rmsdev,iscor
1211 #endif
1212         enddo ! i
1213         write (iout,*) iii," conformations (from",indstart(j)," to",&
1214          indend(j),") read from ",&
1215          bxname(:ilen(bxname))
1216         close (ientin,status="delete")
1217       enddo ! j
1218       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
1219 #if (defined(AIX) && !defined(JUBL))
1220       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
1221 #else
1222       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
1223 #endif
1224 #endif
1225       return
1226   101 write (iout,*) "Error in scratchfile."
1227       call flush(iout)
1228       return 1
1229       end subroutine write_dbase
1230 !-------------------------------------------------------------------------------
1231       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
1232 !      implicit none
1233 !      include "DIMENSIONS"
1234 !      include "DIMENSIONS.ZSCOPT"
1235 !      include "DIMENSIONS.FREE"
1236 !      include "DIMENSIONS.COMPAR"
1237 #ifdef MPI
1238       include "mpif.h"
1239       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1240 !      include "COMMON.MPI"
1241 #endif
1242 !      include "COMMON.CONTROL"
1243 !      include "COMMON.CHAIN"
1244 !      include "COMMON.IOUNITS"
1245 !      include "COMMON.PROTFILES"
1246 !      include "COMMON.NAMES"
1247 !      include "COMMON.VAR"
1248 !      include "COMMON.SBRIDGE"
1249 !      include "COMMON.GEO"
1250 !      include "COMMON.FFIELD"
1251 !      include "COMMON.ENEPS"
1252 !      include "COMMON.LOCAL"
1253 !      include "COMMON.WEIGHTS"
1254 !      include "COMMON.INTERACT"
1255 !      include "COMMON.FREE"
1256 !      include "COMMON.ENERGIES"
1257 !      include "COMMON.COMPAR"
1258 !      include "COMMON.PROT"
1259       integer :: i,j,itmp,iscor,iret,ixdrf
1260       real(kind=8) :: rmsdev,efree,eini
1261       real(kind=4) :: csingle(3,nres*2),xoord(3,2*nres+2)
1262       real(kind=4) :: prec
1263
1264 !      write (iout,*) "cxwrite"
1265 !      call flush(iout)
1266       prec=10000.0
1267       do i=1,nres
1268        do j=1,3
1269         xoord(j,i)=csingle(j,i)
1270        enddo
1271       enddo
1272       do i=nnt,nct
1273        do j=1,3
1274         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
1275        enddo
1276       enddo
1277
1278       itmp=nres+nct-nnt+1
1279
1280 !      write (iout,*) "itmp",itmp
1281 !      call flush(iout)
1282 #if (defined(AIX) && !defined(JUBL))
1283       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
1284
1285 !      write (iout,*) "xdrf3dfcoord"
1286 !      call flush(iout)
1287       call xdrfint_(ixdrf, nss, iret)
1288             write (iout,*) "iret",iret
1289             write (iout,*) "nss",nss,i,"TUTU"
1290       do j=1,nss
1291         call xdrfint_(ixdrf, ihpb(j), iret)
1292         call xdrfint_(ixdrf, jhpb(j), iret)
1293             write(iout,*), ihpb(j),jhpb(j),"TUTU"
1294       enddo
1295       call xdrffloat_(ixdrf,real(eini),iret)
1296       call xdrffloat_(ixdrf,real(efree),iret)
1297             write(iout,*) "TUTU", eini
1298             write(iout,*) "TUTU", efree
1299       call xdrffloat_(ixdrf,real(rmsdev),iret)
1300       call xdrfint_(ixdrf,iscor,iret)
1301 #else
1302       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
1303             write (iout,*) "iret",iret
1304             write (iout,*) "nss",nss,i,"TUTU"
1305
1306       call xdrfint(ixdrf, nss, iret)
1307       do j=1,nss
1308         call xdrfint(ixdrf, ihpb(j), iret)
1309         call xdrfint(ixdrf, jhpb(j), iret)
1310             write(iout,*), ihpb(j),jhpb(j),"TUTU"
1311       enddo
1312       call xdrffloat(ixdrf,real(eini),iret)
1313       call xdrffloat(ixdrf,real(efree),iret)
1314             write(iout,*) "TUTU", eini
1315             write(iout,*) "TUTU", efree
1316       call xdrffloat(ixdrf,real(rmsdev),iret)
1317       call xdrfint(ixdrf,iscor,iret)
1318 #endif
1319
1320       return
1321       end subroutine cxwrite
1322 !-------------------------------------------------------------------------------
1323 ! slices.F
1324 !-------------------------------------------------------------------------------
1325       subroutine set_slices(is,ie,ts,te,iR,ib,iparm)
1326 !      implicit none
1327 !      include 'DIMENSIONS'
1328 !      include 'DIMENSIONS.ZSCOPT'
1329 !      include 'DIMENSIONS.FREE'
1330 !      include 'COMMON.IOUNITS'
1331 !      include 'COMMON.PROTFILES'
1332 !      include 'COMMON.OBCINKA'
1333 !      include 'COMMON.PROT'
1334       integer :: islice,iR,ib,iparm
1335       integer :: is(MaxSlice),ie(MaxSlice),nrec_slice
1336       real(kind=8) :: ts(MaxSlice),te(MaxSlice),time_slice
1337       time_slice=0
1338       do islice=1,nslice
1339         if (time_end_collect(iR,ib,iparm).ge.1.0d10) then
1340           ts(islice)=time_start_collect(iR,ib,iparm)
1341           te(islice)=time_end_collect(iR,ib,iparm)
1342           nrec_slice=(rec_end(iR,ib,iparm)- &
1343              rec_start(iR,ib,iparm)+1)/nslice
1344           is(islice)=rec_start(iR,ib,iparm)+(islice-1)*nrec_slice
1345           ie(islice)=rec_start(iR,ib,iparm)+islice*nrec_slice-1
1346         else
1347           time_slice=(time_end_collect(iR,ib,iparm) &
1348           -time_start_collect(iR,ib,iparm))/nslice
1349           ts(islice)=time_start_collect(iR,ib,iparm)+(islice-1)* &
1350            time_slice
1351           te(islice)=time_start_collect(iR,ib,iparm)+islice*time_slice
1352           is(islice)=rec_start(iR,ib,iparm)
1353           ie(islice)=rec_end(iR,ib,iparm)
1354         endif
1355       enddo
1356
1357       write (iout,*) "nrec_slice",nrec_slice," time_slice",time_slice
1358       write (iout,*) "is",(is(islice),islice=1,nslice)
1359       write (iout,*) "ie",(ie(islice),islice=1,nslice)
1360       write (iout,*) "rec_start",&
1361         rec_start(iR,ib,iparm)," rec_end",rec_end(iR,ib,iparm)
1362       write (iout,*) "ts",(ts(islice),islice=1,nslice)
1363       write (iout,*) "te",(te(islice),islice=1,nslice)
1364       write (iout,*) "time_start",&
1365         time_start_collect(iR,ib,iparm)," time_end",&
1366         time_end_collect(iR,ib,iparm)
1367       call flush(iout)
1368
1369       return
1370       end subroutine set_slices
1371 !-----------------------------------------------------------------------------
1372       integer function slice(irecord,time,is,ie,ts,te)
1373 !      implicit none
1374 !      include 'DIMENSIONS'
1375 !      include 'DIMENSIONS.ZSCOPT'
1376 !      include 'DIMENSIONS.FREE'
1377 !      include 'COMMON.IOUNITS'
1378 !      include 'COMMON.PROTFILES'
1379 !      include 'COMMON.OBCINKA'
1380 !      include 'COMMON.PROT'
1381       integer :: is(MaxSlice),ie(MaxSlice),nrec_slice
1382       real(kind=8) :: ts(MaxSlice),te(MaxSlice),time_slice
1383       integer :: i,ii,irecord
1384       real(kind=8) :: time
1385
1386 !      write (iout,*) "within slice nslice",nslice
1387 !      call flush(iout)
1388       if (irecord.lt.is(1) .or. time.lt.ts(1)) then
1389         ii=0
1390       else
1391         ii=1
1392         do while (ii.le.nslice .and. &
1393                  (irecord.lt.is(ii) .or. irecord.gt.ie(ii) .or. &
1394                  time.lt.ts(ii) .or. time.gt.te(ii)) )
1395 !          write (iout,*) "ii",ii,time,ts(ii)
1396 !          call flush(iout)
1397           ii=ii+1
1398         enddo
1399       endif
1400 !      write (iout,*) "end: ii",ii
1401 !      call flush(iout)
1402       slice=ii
1403       return
1404       end function slice
1405 !-----------------------------------------------------------------------------
1406 ! enecalc1.F
1407 !-----------------------------------------------------------------------------
1408       logical function conf_check(ii,iprint)
1409
1410       use names, only:ntyp1
1411       use geometry_data
1412       use energy_data, only:itype,dsc,molnum
1413       use geometry, only:int_from_cart1
1414 !      use
1415 !      include "DIMENSIONS"
1416 !      include "DIMENSIONS.ZSCOPT"
1417 !      include "DIMENSIONS.FREE"
1418 !#ifdef MPI
1419 !      use MPI_data
1420 !      include "mpif.h"
1421 !      include "COMMON.MPI"
1422 !#endif
1423 !      include "COMMON.CHAIN"
1424 !      include "COMMON.IOUNITS"
1425 !      include "COMMON.PROTFILES"
1426 !      include "COMMON.NAMES"
1427 !      include "COMMON.VAR"
1428 !      include "COMMON.SBRIDGE"
1429 !      include "COMMON.GEO"
1430 !      include "COMMON.FFIELD"
1431 !      include "COMMON.ENEPS"
1432 !      include "COMMON.LOCAL"
1433 !      include "COMMON.WEIGHTS"
1434 !      include "COMMON.INTERACT"
1435 !      include "COMMON.FREE"
1436 !      include "COMMON.ENERGIES"
1437 !      include "COMMON.CONTROL"
1438 !      include "COMMON.TORCNSTR"
1439 !      implicit none
1440 #ifdef MPI
1441       include "mpif.h"
1442       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1443 #endif
1444       integer :: j,k,l,ii,itj,iprint,mnum
1445       if (.not. check_conf) then
1446         conf_check=.true.
1447         return
1448       endif
1449       call int_from_cart1(.false.)
1450       do j=nnt+1,nct
1451          mnum=molnum(j)
1452          if (mnum.ne.1) cycle
1453          if (mnum.eq.5) cycle
1454         if (itype(j-1,mnum).ne.ntyp1 .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
1455           (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
1456           if (iprint.gt.0) &
1457           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
1458             " for conformation",ii,mnum
1459           if (iprint.gt.1) then
1460             write (iout,*) "The Cartesian geometry is:"
1461             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1462             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1463             write (iout,*) "The internal geometry is:"
1464             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1465             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1466             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1467             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1468             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1469             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1470           endif
1471           if (iprint.gt.0) write (iout,*) &
1472             "This conformation WILL NOT be added to the database."
1473           conf_check=.false.
1474           return
1475         endif
1476       enddo
1477       do j=nnt,nct
1478         mnum=molnum(j)
1479         if (mnum.ne.1) cycle
1480         itj=itype(j,mnum)
1481         if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
1482            (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
1483           if (iprint.gt.0) &
1484           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
1485            " for conformation",ii
1486           if (iprint.gt.1) then
1487             write (iout,*) "The Cartesian geometry is:"
1488             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1489             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1490             write (iout,*) "The internal geometry is:"
1491             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1492             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1493             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1494             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1495             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1496             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1497           endif
1498           if (iprint.gt.0) write (iout,*) &
1499             "This conformation WILL NOT be added to the database."
1500           conf_check=.false.
1501           return
1502         endif
1503       enddo
1504       do j=3,nres
1505         mnum=molnum(j)
1506         itj=itype(j,mnum)
1507         if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
1508         if (itype(j-1,mnum).eq.ntyp1_molec(mnum)) cycle
1509         if (itype(j-2,mnum).eq.ntyp1_molec(mnum)) cycle
1510         if (theta(j).le.0.0d0) then
1511           if (iprint.gt.0) &
1512           write (iout,*) "Zero theta angle(s) in conformation",ii
1513           if (iprint.gt.1) then
1514             write (iout,*) "The Cartesian geometry is:"
1515             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1516             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1517             write (iout,*) "The internal geometry is:"
1518             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1519             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1520             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1521             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1522             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1523             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1524           endif
1525           if (iprint.gt.0) write (iout,*) &
1526             "This conformation WILL NOT be added to the database."
1527           conf_check=.false.
1528           return
1529         endif
1530         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
1531       enddo
1532       conf_check=.true.
1533 !      write (iout,*) "conf_check passed",ii
1534       return
1535       end function conf_check
1536 !-----------------------------------------------------------------------------
1537       end module io_database