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
 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
 920 !      implicit none
 921 !      include "DIMENSIONS"
 922 !      include "DIMENSIONS.ZSCOPT"
 923 !      include "DIMENSIONS.FREE"
 924 !      include "DIMENSIONS.COMPAR"
 925       use geometry, only:int_from_cart1
 926 #ifdef MPI
 927       include "mpif.h"
 928       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 929 !      include "COMMON.MPI"
 930 #endif
 931 !      include "COMMON.CONTROL"
 932 !      include "COMMON.CHAIN"
 933 !      include "COMMON.IOUNITS"
 934 !      include "COMMON.PROTFILES"
 935 !      include "COMMON.NAMES"
 936 !      include "COMMON.VAR"
 937 !      include "COMMON.SBRIDGE"
 938 !      include "COMMON.GEO"
 939 !      include "COMMON.FFIELD"
 940 !      include "COMMON.ENEPS"
 941 !      include "COMMON.LOCAL"
 942 !      include "COMMON.WEIGHTS"
 943 !      include "COMMON.INTERACT"
 944 !      include "COMMON.FREE"
 945 !      include "COMMON.ENERGIES"
 946 !      include "COMMON.COMPAR"
 947 !      include "COMMON.PROT"
 948 !      include "COMMON.CONTACTS1"
 949       character(len=64) :: nazwa
 950       character(len=80) :: bxname,cxname
 951       character(len=64) :: bprotfile_temp
 952       character(len=3) :: liczba,licz
 953       character(len=2) :: licz2
 954       integer :: i,itj,ii,iii,j,k,l
 955       integer :: ixdrf,iret
 956       integer :: iscor,islice
 957       real(kind=8) :: rmsdev,efree,eini
 958       real(kind=4) :: csingle(3,nres*2)
 959       real(kind=8) :: energ
 960 !      integer ilen,iroof
 961 !      external ilen,iroof
 962       integer :: ir,ib,iparm
 963       integer :: isecstr(nres)
 964       write (licz2,'(bz,i2.2)') islice
 965       call opentmp(islice,ientout,bprotfile_temp)
 966       write (iout,*) "bprotfile_temp ",bprotfile_temp
 967       call flush(iout)
 968       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0 &
 969          .and. ensembles.eq.0) then
 970         close(ientout,status="delete")
 971         return
 972       endif
 973 #ifdef MPI
 974       write (liczba,'(bz,i3.3)') me
 975       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 976         if (.not.separate_parset) then
 977           bxname = prefix(:ilen(prefix))//liczba//".bx"
 978         else
 979           write (licz,'(bz,i3.3)') myparm
 980           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 981         endif
 982         print *,bxname
 983         open (ientin,file=bxname,status="unknown",&
 984           form="unformatted",access="direct",recl=lenrec1)
 985       endif
 986 #else
 987       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 988         if (nslice.eq.1) then
 989           bxname = prefix(:ilen(prefix))//".bx"
 990         else
 991           bxname = prefix(:ilen(prefix))// &
 992                  "_slice_"//licz2//".bx"
 993         endif
 994         open (ientin,file=bxname,status="unknown",&
 995           form="unformatted",access="direct",recl=lenrec1)
 996         write (iout,*) "Calculating energies; writing geometry",&
 997        " and energy components to ",bxname(:ilen(bxname))
 998       endif
 999 #if (defined(AIX) && !defined(JUBL))
1000         call xdrfopen_(ixdrf,cxname, "w", iret)
1001 #else
1002         call xdrfopen(ixdrf,cxname, "w", iret)
1003 #endif
1004         if (iret.eq.0) then
1005           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
1006           cxfile=.false.
1007         endif
1008 !el      endif
1009 #endif
1010       print *,indpdb
1011       if (indpdb.gt.0) then
1012         if (nslice.eq.1) then
1013 #ifdef MPI
1014          if (.not.separate_parset) then
1015            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot)) &
1016              //liczba//'.stat'
1017          else
1018            write (licz,'(bz,i3.3)') myparm
1019            statname=prefix(:ilen(prefix))//'_par'//licz//'_'// &
1020             pot(:ilen(pot))//liczba//'.stat'
1021          endif
1022          print *,statname
1023 #else
1024           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
1025 #endif
1026         else
1027 #ifdef MPI
1028          if (.not.separate_parset) then
1029           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))// &
1030             "_slice_"//licz2//liczba//'.stat'
1031          else
1032           write (licz,'(bz,i3.3)') myparm
1033           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))// &
1034             '_par'//licz//"_slice_"//licz2//liczba//'.stat'
1035          endif
1036 #else
1037           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot)) &
1038             //"_slice_"//licz2//'.stat'
1039 #endif
1040         endif
1041         print *,istat,statname
1042         open(istat,file=statname,status="unknown")
1043       endif
1044
1045 #ifdef MPI
1046       do i=1,scount(me)
1047 #else
1048       do i=1,ntot(islice)
1049 #endif
1050         read(ientout,rec=i,err=101) &
1051           ((csingle(l,k),l=1,3),k=1,nres),&
1052           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1053           nss,(ihpb(k),jhpb(k),k=1,nss),&
1054           eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
1055 !        write (iout,*) iR,ib,iparm,eini,efree
1056         do j=1,2*nres
1057           do k=1,3
1058             c(k,j)=csingle(k,j)
1059           enddo
1060         enddo
1061         call int_from_cart1(.false.)
1062         iscore=0
1063 !        write (iout,*) "Calling conf_compar",i
1064 !        call flush(iout)
1065          anatemp= 1.0d0/(beta_h(ib,iparm)*1.987D-3)
1066         if (indpdb.gt.0) then
1067           call conf_compar(i,.false.,.true.)
1068 !        else
1069 !            call elecont(.false.,ncont,icont,nnt,nct)
1070 !            call secondary2(.false.,.false.,ncont,icont,isecstr)
1071         endif
1072 !        write (iout,*) "Exit conf_compar",i
1073 !        call flush(iout)
1074         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i) &
1075           ((csingle(l,k),l=1,3),k=1,nres),&
1076           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1077           nss,(ihpb(k),jhpb(k),k=1,nss),&
1078 !     &    potE(i,iparm),-entfac(i),rms_nat,iscore
1079           potE(i,nparmset),-entfac(i),rms_nat,iscore
1080 !        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
1081 #ifndef MPI
1082         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),&
1083           -entfac(i),rms_nat,iscore)
1084 #endif
1085       enddo
1086       close(ientout,status="delete")
1087       close(istat)
1088       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
1089 #ifdef MPI
1090       call MPI_Barrier(WHAM_COMM,IERROR)
1091       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile &
1092          .and. ensembles.eq.0) return
1093       write (iout,*)
1094       if (bxfile .or. ensembles.gt.0) then
1095         if (nslice.eq.1) then
1096           if (.not.separate_parset) then
1097             bxname = prefix(:ilen(prefix))//".bx"
1098           else
1099             write (licz,'(bz,i3.3)') myparm
1100             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
1101           endif
1102         else
1103           if (.not.separate_parset) then
1104             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
1105           else
1106             write (licz,'(bz,i3.3)') myparm
1107             bxname = prefix(:ilen(prefix))//"par_"//licz// &
1108               "_slice_"//licz2//".bx"
1109           endif
1110         endif
1111         open (ientout,file=bxname,status="unknown",&
1112             form="unformatted",access="direct",recl=lenrec1)
1113         write (iout,*) "Master is creating binary database ",&
1114          bxname(:ilen(bxname))
1115       endif
1116       if (cxfile) then
1117         if (nslice.eq.1) then
1118           if (.not.separate_parset) then
1119             cxname = prefix(:ilen(prefix))//".cx"
1120           else
1121             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
1122           endif
1123         else
1124           if (.not.separate_parset) then
1125             cxname = prefix(:ilen(prefix))// &
1126                    "_slice_"//licz2//".cx"
1127           else
1128             cxname = prefix(:ilen(prefix))//"_par"//licz// &
1129                    "_slice_"//licz2//".cx"
1130           endif
1131         endif
1132 #if (defined(AIX) && !defined(JUBL))
1133         call xdrfopen_(ixdrf,cxname, "w", iret)
1134 #else
1135         call xdrfopen(ixdrf,cxname, "w", iret)
1136 #endif
1137         if (iret.eq.0) then
1138           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
1139           cxfile=.false.
1140         endif
1141       endif
1142       do j=0,nprocs-1
1143         write (liczba,'(bz,i3.3)') j
1144         if (separate_parset) then
1145           write (licz,'(bz,i3.3)') myparm
1146           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
1147         else
1148           bxname = prefix(:ilen(prefix))//liczba//".bx"
1149         endif
1150         open (ientin,file=bxname,status="unknown",&
1151           form="unformatted",access="direct",recl=lenrec1)
1152         write (iout,*) "Master is reading conformations from ",&
1153          bxname(:ilen(bxname))
1154         iii = 0
1155 !        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
1156 !        call flush(iout)
1157         do i=indstart(j),indend(j)
1158           iii = iii+1
1159           read(ientin,rec=iii,err=101) &
1160             ((csingle(l,k),l=1,3),k=1,nres),&
1161             ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1162             nss,(ihpb(k),jhpb(k),k=1,nss),&
1163             eini,efree,rmsdev,iscor
1164           if (bxfile .or. ensembles.gt.0) then
1165             write (ientout,rec=i) &
1166               ((csingle(l,k),l=1,3),k=1,nres),&
1167               ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
1168               nss,(ihpb(k),jhpb(k),k=1,nss),&
1169               eini,efree,rmsdev,iscor
1170           endif
1171            print *,"before cxwrite"
1172           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
1173           print *,"after cxwrite"
1174 #ifdef DEBUG
1175           do k=1,2*nres
1176             do l=1,3
1177               c(l,k)=csingle(l,k)
1178             enddo
1179           enddo
1180           call int_from_cart1(.false.)
1181           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
1182           write (iout,*) "The Cartesian geometry is:"
1183           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1184           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1185           write (iout,*) "The internal geometry is:"
1186           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1187           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1188           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1189           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1190           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1191           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1192           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
1193           write (iout,'(f10.5,i5)') rmsdev,iscor
1194 #endif
1195         enddo ! i
1196         write (iout,*) iii," conformations (from",indstart(j)," to",&
1197          indend(j),") read from ",&
1198          bxname(:ilen(bxname))
1199         close (ientin,status="delete")
1200       enddo ! j
1201       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
1202 #if (defined(AIX) && !defined(JUBL))
1203       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
1204 #else
1205       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
1206 #endif
1207 #endif
1208       return
1209   101 write (iout,*) "Error in scratchfile."
1210       call flush(iout)
1211       return 1
1212       end subroutine write_dbase
1213 !-------------------------------------------------------------------------------
1214       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
1215 !      implicit none
1216 !      include "DIMENSIONS"
1217 !      include "DIMENSIONS.ZSCOPT"
1218 !      include "DIMENSIONS.FREE"
1219 !      include "DIMENSIONS.COMPAR"
1220 #ifdef MPI
1221       include "mpif.h"
1222       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1223 !      include "COMMON.MPI"
1224 #endif
1225 !      include "COMMON.CONTROL"
1226 !      include "COMMON.CHAIN"
1227 !      include "COMMON.IOUNITS"
1228 !      include "COMMON.PROTFILES"
1229 !      include "COMMON.NAMES"
1230 !      include "COMMON.VAR"
1231 !      include "COMMON.SBRIDGE"
1232 !      include "COMMON.GEO"
1233 !      include "COMMON.FFIELD"
1234 !      include "COMMON.ENEPS"
1235 !      include "COMMON.LOCAL"
1236 !      include "COMMON.WEIGHTS"
1237 !      include "COMMON.INTERACT"
1238 !      include "COMMON.FREE"
1239 !      include "COMMON.ENERGIES"
1240 !      include "COMMON.COMPAR"
1241 !      include "COMMON.PROT"
1242       integer :: i,j,itmp,iscor,iret,ixdrf
1243       real(kind=8) :: rmsdev,efree,eini
1244       real(kind=4) :: csingle(3,nres*2),xoord(3,2*nres+2)
1245       real(kind=4) :: prec
1246
1247 !      write (iout,*) "cxwrite"
1248 !      call flush(iout)
1249       prec=10000.0
1250       do i=1,nres
1251        do j=1,3
1252         xoord(j,i)=csingle(j,i)
1253        enddo
1254       enddo
1255       do i=nnt,nct
1256        do j=1,3
1257         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
1258        enddo
1259       enddo
1260
1261       itmp=nres+nct-nnt+1
1262
1263 !      write (iout,*) "itmp",itmp
1264 !      call flush(iout)
1265 #if (defined(AIX) && !defined(JUBL))
1266       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
1267
1268 !      write (iout,*) "xdrf3dfcoord"
1269 !      call flush(iout)
1270       call xdrfint_(ixdrf, nss, iret)
1271       do j=1,nss
1272         call xdrfint_(ixdrf, ihpb(j), iret)
1273         call xdrfint_(ixdrf, jhpb(j), iret)
1274       enddo
1275       call xdrffloat_(ixdrf,real(eini),iret)
1276       call xdrffloat_(ixdrf,real(efree),iret)
1277       call xdrffloat_(ixdrf,real(rmsdev),iret)
1278       call xdrfint_(ixdrf,iscor,iret)
1279 #else
1280       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
1281
1282       call xdrfint(ixdrf, nss, iret)
1283       do j=1,nss
1284         call xdrfint(ixdrf, ihpb(j), iret)
1285         call xdrfint(ixdrf, jhpb(j), iret)
1286       enddo
1287       call xdrffloat(ixdrf,real(eini),iret)
1288       call xdrffloat(ixdrf,real(efree),iret)
1289       call xdrffloat(ixdrf,real(rmsdev),iret)
1290       call xdrfint(ixdrf,iscor,iret)
1291 #endif
1292
1293       return
1294       end subroutine cxwrite
1295 !-------------------------------------------------------------------------------
1296 ! slices.F
1297 !-------------------------------------------------------------------------------
1298       subroutine set_slices(is,ie,ts,te,iR,ib,iparm)
1299 !      implicit none
1300 !      include 'DIMENSIONS'
1301 !      include 'DIMENSIONS.ZSCOPT'
1302 !      include 'DIMENSIONS.FREE'
1303 !      include 'COMMON.IOUNITS'
1304 !      include 'COMMON.PROTFILES'
1305 !      include 'COMMON.OBCINKA'
1306 !      include 'COMMON.PROT'
1307       integer :: islice,iR,ib,iparm
1308       integer :: is(MaxSlice),ie(MaxSlice),nrec_slice
1309       real(kind=8) :: ts(MaxSlice),te(MaxSlice),time_slice
1310
1311       do islice=1,nslice
1312         if (time_end_collect(iR,ib,iparm).ge.1.0d10) then
1313           ts(islice)=time_start_collect(iR,ib,iparm)
1314           te(islice)=time_end_collect(iR,ib,iparm)
1315           nrec_slice=(rec_end(iR,ib,iparm)- &
1316              rec_start(iR,ib,iparm)+1)/nslice
1317           is(islice)=rec_start(iR,ib,iparm)+(islice-1)*nrec_slice
1318           ie(islice)=rec_start(iR,ib,iparm)+islice*nrec_slice-1
1319         else
1320           time_slice=(time_end_collect(iR,ib,iparm) &
1321           -time_start_collect(iR,ib,iparm))/nslice
1322           ts(islice)=time_start_collect(iR,ib,iparm)+(islice-1)* &
1323            time_slice
1324           te(islice)=time_start_collect(iR,ib,iparm)+islice*time_slice
1325           is(islice)=rec_start(iR,ib,iparm)
1326           ie(islice)=rec_end(iR,ib,iparm)
1327         endif
1328       enddo
1329
1330       write (iout,*) "nrec_slice",nrec_slice," time_slice",time_slice
1331       write (iout,*) "is",(is(islice),islice=1,nslice)
1332       write (iout,*) "ie",(ie(islice),islice=1,nslice)
1333       write (iout,*) "rec_start",&
1334         rec_start(iR,ib,iparm)," rec_end",rec_end(iR,ib,iparm)
1335       write (iout,*) "ts",(ts(islice),islice=1,nslice)
1336       write (iout,*) "te",(te(islice),islice=1,nslice)
1337       write (iout,*) "time_start",&
1338         time_start_collect(iR,ib,iparm)," time_end",&
1339         time_end_collect(iR,ib,iparm)
1340       call flush(iout)
1341
1342       return
1343       end subroutine set_slices
1344 !-----------------------------------------------------------------------------
1345       integer function slice(irecord,time,is,ie,ts,te)
1346 !      implicit none
1347 !      include 'DIMENSIONS'
1348 !      include 'DIMENSIONS.ZSCOPT'
1349 !      include 'DIMENSIONS.FREE'
1350 !      include 'COMMON.IOUNITS'
1351 !      include 'COMMON.PROTFILES'
1352 !      include 'COMMON.OBCINKA'
1353 !      include 'COMMON.PROT'
1354       integer :: is(MaxSlice),ie(MaxSlice),nrec_slice
1355       real(kind=8) :: ts(MaxSlice),te(MaxSlice),time_slice
1356       integer :: i,ii,irecord
1357       real(kind=8) :: time
1358
1359 !      write (iout,*) "within slice nslice",nslice
1360 !      call flush(iout)
1361       if (irecord.lt.is(1) .or. time.lt.ts(1)) then
1362         ii=0
1363       else
1364         ii=1
1365         do while (ii.le.nslice .and. &
1366                  (irecord.lt.is(ii) .or. irecord.gt.ie(ii) .or. &
1367                  time.lt.ts(ii) .or. time.gt.te(ii)) )
1368 !          write (iout,*) "ii",ii,time,ts(ii)
1369 !          call flush(iout)
1370           ii=ii+1
1371         enddo
1372       endif
1373 !      write (iout,*) "end: ii",ii
1374 !      call flush(iout)
1375       slice=ii
1376       return
1377       end function slice
1378 !-----------------------------------------------------------------------------
1379 ! enecalc1.F
1380 !-----------------------------------------------------------------------------
1381       logical function conf_check(ii,iprint)
1382
1383       use names, only:ntyp1
1384       use geometry_data
1385       use energy_data, only:itype,dsc,molnum
1386       use geometry, only:int_from_cart1
1387 !      use
1388 !      include "DIMENSIONS"
1389 !      include "DIMENSIONS.ZSCOPT"
1390 !      include "DIMENSIONS.FREE"
1391 !#ifdef MPI
1392 !      use MPI_data
1393 !      include "mpif.h"
1394 !      include "COMMON.MPI"
1395 !#endif
1396 !      include "COMMON.CHAIN"
1397 !      include "COMMON.IOUNITS"
1398 !      include "COMMON.PROTFILES"
1399 !      include "COMMON.NAMES"
1400 !      include "COMMON.VAR"
1401 !      include "COMMON.SBRIDGE"
1402 !      include "COMMON.GEO"
1403 !      include "COMMON.FFIELD"
1404 !      include "COMMON.ENEPS"
1405 !      include "COMMON.LOCAL"
1406 !      include "COMMON.WEIGHTS"
1407 !      include "COMMON.INTERACT"
1408 !      include "COMMON.FREE"
1409 !      include "COMMON.ENERGIES"
1410 !      include "COMMON.CONTROL"
1411 !      include "COMMON.TORCNSTR"
1412 !      implicit none
1413 #ifdef MPI
1414       include "mpif.h"
1415       integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1416 #endif
1417       integer :: j,k,l,ii,itj,iprint,mnum
1418       if (.not. check_conf) then
1419         conf_check=.true.
1420         return
1421       endif
1422       call int_from_cart1(.false.)
1423       do j=nnt+1,nct
1424          mnum=molnum(j)
1425          if (mnum.ne.1) cycle
1426          if (mnum.eq.5) cycle
1427         if (itype(j-1,mnum).ne.ntyp1 .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
1428           (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
1429           if (iprint.gt.0) &
1430           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),&
1431             " for conformation",ii,mnum
1432           if (iprint.gt.1) then
1433             write (iout,*) "The Cartesian geometry is:"
1434             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1435             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1436             write (iout,*) "The internal geometry is:"
1437             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1438             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1439             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1440             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1441             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1442             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1443           endif
1444           if (iprint.gt.0) write (iout,*) &
1445             "This conformation WILL NOT be added to the database."
1446           conf_check=.false.
1447           return
1448         endif
1449       enddo
1450       do j=nnt,nct
1451         mnum=molnum(j)
1452         if (mnum.ne.1) cycle
1453         itj=itype(j,mnum)
1454         if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1_molec(mnum) .and. &
1455            (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
1456           if (iprint.gt.0) &
1457           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),&
1458            " for conformation",ii
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=3,nres
1478         if (theta(j).le.0.0d0) then
1479           if (iprint.gt.0) &
1480           write (iout,*) "Zero theta angle(s) in conformation",ii
1481           if (iprint.gt.1) then
1482             write (iout,*) "The Cartesian geometry is:"
1483             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1484             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1485             write (iout,*) "The internal geometry is:"
1486             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1487             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1488             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1489             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1490             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1491             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1492           endif
1493           if (iprint.gt.0) write (iout,*) &
1494             "This conformation WILL NOT be added to the database."
1495           conf_check=.false.
1496           return
1497         endif
1498         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
1499       enddo
1500       conf_check=.true.
1501 !      write (iout,*) "conf_check passed",ii
1502       return
1503       end function conf_check
1504 !-----------------------------------------------------------------------------
1505       end module io_database