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