source/wham/src-M/readrtns.F

   1       subroutine read_general_data(*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       include "DIMENSIONS.FREE"
   6       include "COMMON.TORSION"
   7       include "COMMON.INTERACT"
   8       include "COMMON.IOUNITS"
   9       include "COMMON.TIME1"
  10       include "COMMON.PROT"
  11       include "COMMON.PROTFILES"
  12       include "COMMON.CHAIN"
  13       include "COMMON.NAMES"
  14       include "COMMON.FFIELD"
  15       include "COMMON.ENEPS"
  16       include "COMMON.WEIGHTS"
  17       include "COMMON.FREE"
  18       include "COMMON.CONTROL"
  19       include "COMMON.ENERGIES"
  20       include "COMMON.SPLITELE"
  21       include "COMMON.SBRIDGE"
  22       include "COMMON.SHIELD"
  23       character*800 controlcard
  24       integer i,j,k,ii,n_ene_found
  25       integer ind,itype1,itype2,itypf,itypsc,itypp
  26       integer ilen
  27       external ilen
  28       character*16 ucase
  29       character*16 key
  30       external ucase
  31       double precision pi
  32       call card_concat(controlcard,.true.)
  33       call readi(controlcard,"N_ENE",n_ene,max_ene)
  34       if (n_ene.gt.max_ene) then
  35         write (iout,*) "Error: parameter out of range: N_ENE",n_ene,
  36      &    max_ene
  37         return1
  38       endif
  39       call readi(controlcard,"NPARMSET",nparmset,1)
  40       separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
  41       call readi(controlcard,"IPARMPRINT",iparmprint,1)
  42       write (iout,*) "PARMPRINT",iparmprint
  43       if (nparmset.gt.max_parm) then
  44         write (iout,*) "Error: parameter out of range: NPARMSET",
  45      &    nparmset, Max_Parm
  46         return1
  47       endif
  48       call readi(controlcard,"MAXIT",maxit,5000)
  49       call reada(controlcard,"FIMIN",fimin,1.0d-3)
  50       call readi(controlcard,"ENSEMBLES",ensembles,0)
  51       hamil_rep=index(controlcard,"HAMIL_REP").gt.0
  52       write (iout,*) "Number of energy parameter sets",nparmset
  53       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
  54       write (iout,*) "MaxSlice",MaxSlice
  55       call readi(controlcard,"NSLICE",nslice,1)
  56       call flush(iout)
  57       if (nslice.gt.MaxSlice) then
  58         write (iout,*) "Error: parameter out of range: NSLICE",nslice,
  59      &    MaxSlice
  60         return1
  61       endif
  62       write (iout,*) "Frequency of storing conformations",
  63      & (isampl(i),i=1,nparmset)
  64       write (iout,*) "Maxit",maxit," Fimin",fimin
  65       call readi(controlcard,"NQ",nQ,1)
  66       if (nQ.gt.MaxQ) then
  67         write (iout,*) "Error: parameter out of range: NQ",nq,
  68      &    maxq
  69         return1
  70       endif
  71       indpdb=0
  72       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
  73       call reada(controlcard,"DELTA",delta,1.0d-2)
  74       call readi(controlcard,"EINICHECK",einicheck,2)
  75       call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
  76       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
  77       call readi(controlcard,"RESCALE",rescale_mode,1)
  78       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
  79       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
  80        call reada(controlcard,'BOXX',boxxsize,100.0d0)
  81        call reada(controlcard,'BOXY',boxysize,100.0d0)
  82        call reada(controlcard,'BOXZ',boxzsize,100.0d0)
  83       call readi(controlcard,'TUBEMOD',tubelog,0)
  84       write (iout,*) TUBElog,"TUBEMODE"
  85       call readi(controlcard,'GENCONSTR',genconstr,0)
  86
  87 c Cutoff range for interactions
  88       call reada(controlcard,"R_CUT",r_cut,15.0d0)
  89       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
  90       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
  91       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
  92       if (lipthick.gt.0.0d0) then
  93        bordliptop=(boxzsize+lipthick)/2.0
  94        bordlipbot=bordliptop-lipthick
  95 C      endif
  96       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
  97      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
  98       buflipbot=bordlipbot+lipbufthick
  99       bufliptop=bordliptop-lipbufthick
 100       if ((lipbufthick*2.0d0).gt.lipthick)
 101      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 102       endif
 103       write(iout,*) "bordliptop=",bordliptop
 104       write(iout,*) "bordlipbot=",bordlipbot
 105       write(iout,*) "bufliptop=",bufliptop
 106       write(iout,*) "buflipbot=",buflipbot
 107       call readi(controlcard,'SYM',symetr,1)
 108       write (iout,*) "DISTCHAINMAX",distchainmax
 109       write (iout,*) "delta",delta
 110       write (iout,*) "einicheck",einicheck
 111       write (iout,*) "rescale_mode",rescale_mode
 112       call flush(iout)
 113       if (TUBElog.gt.0) then
 114        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 115        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 116        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
 117        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 118        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 119        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 120        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 121        buftubebot=bordtubebot+tubebufthick
 122        buftubetop=bordtubetop-tubebufthick
 123       endif
 124       bxfile=index(controlcard,"BXFILE").gt.0
 125       cxfile=index(controlcard,"CXFILE").gt.0
 126       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
 127      & bxfile=.true.
 128       histfile=index(controlcard,"HISTFILE").gt.0
 129       histout=index(controlcard,"HISTOUT").gt.0
 130       entfile=index(controlcard,"ENTFILE").gt.0
 131       zscfile=index(controlcard,"ZSCFILE").gt.0
 132       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
 133       write (iout,*) "with_dihed_constr ",with_dihed_constr
 134       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 135       write (iout,*) "with_theta_constr ",with_theta_constr
 136       call readi(controlcard,'SHIELD',shield_mode,0)
 137 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 138         write(iout,*) "shield_mode",shield_mode
 139 C      endif
 140       call readi(controlcard,'TORMODE',tor_mode,0)
 141 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 142         write(iout,*) "torsional and valence angle mode",tor_mode
 143       if (shield_mode.gt.0) then
 144       pi=3.141592d0
 145 C VSolvSphere the volume of solving sphere
 146 C      print *,pi,"pi"
 147 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 148 C there will be no distinction between proline peptide group and normal peptide
 149 C group in case of shielding parameters
 150       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 151       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 152       write (iout,*) VSolvSphere,VSolvSphere_div
 153 C long axis of side chain
 154 C      do i=1,ntyp
 155 C      long_r_sidechain(i)=vbldsc0(1,i)
 156 C      short_r_sidechain(i)=sigma0(i)
 157 C      enddo
 158       buff_shield=1.0d0
 159       endif
 160
 161       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 162       dyn_ss=index(controlcard,"DYN_SS").gt.0
 163       return
 164       end
 165 c------------------------------------------------------------------------------
 166       subroutine read_efree(*)
 167 C
 168 C Read molecular data
 169 C
 170       implicit none
 171       include 'DIMENSIONS'
 172       include 'DIMENSIONS.ZSCOPT'
 173       include 'DIMENSIONS.COMPAR'
 174       include 'DIMENSIONS.FREE'
 175       include 'COMMON.IOUNITS'
 176       include 'COMMON.TIME1'
 177       include 'COMMON.SBRIDGE'
 178       include 'COMMON.CONTROL'
 179       include 'COMMON.CHAIN'
 180       include 'COMMON.HEADER'
 181       include 'COMMON.GEO'
 182       include 'COMMON.FREE'
 183       character*320 controlcard,ucase
 184       integer iparm,ib,i,j,npars
 185       integer ilen
 186       external ilen
 187
 188       if (hamil_rep) then
 189         npars=1
 190       else
 191         npars=nParmSet
 192       endif
 193
 194       do iparm=1,npars
 195
 196       call card_concat(controlcard,.true.)
 197       call readi(controlcard,'NT',nT_h(iparm),1)
 198       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
 199       call flush(iout)
 200       if (nT_h(iparm).gt.MaxT_h) then
 201         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),
 202      &    MaxT_h
 203         return1
 204       endif
 205       replica(iparm)=index(controlcard,"REPLICA").gt.0
 206       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
 207       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
 208       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
 209      &  replica(iparm)," umbrella ",umbrella(iparm),
 210      &  " read_iset",read_iset(iparm)
 211       call flush(iout)
 212       do ib=1,nT_h(iparm)
 213         call card_concat(controlcard,.true.)
 214         call readi(controlcard,'NR',nR(ib,iparm),1)
 215         if (umbrella(iparm)) then
 216           nRR(ib,iparm)=1
 217         else
 218           nRR(ib,iparm)=nR(ib,iparm)
 219         endif
 220         if (nR(ib,iparm).gt.MaxR) then
 221           write (iout,*)  "Error: parameter out of range: NR",
 222      &      nR(ib,iparm),MaxR
 223         return1
 224         endif
 225         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
 226         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
 227         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
 228      &    0.0d0)
 229         do i=1,nR(ib,iparm)
 230           call card_concat(controlcard,.true.)
 231           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
 232      &      100.0d0)
 233           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
 234      &      0.0d0)
 235         enddo
 236       enddo
 237       do ib=1,nT_h(iparm)
 238         write (iout,*) "ib",ib," beta_h",
 239      &    1.0d0/(0.001987*beta_h(ib,iparm))
 240         write (iout,*) "nR",nR(ib,iparm)
 241         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
 242         do i=1,nR(ib,iparm)
 243           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
 244      &      "q0",(q0(j,i,ib,iparm),j=1,nQ)
 245         enddo
 246         call flush(iout)
 247       enddo
 248
 249       enddo
 250
 251       if (hamil_rep) then
 252
 253        do iparm=2,nParmSet
 254           nT_h(iparm)=nT_h(1)
 255          do ib=1,nT_h(iparm)
 256            nR(ib,iparm)=nR(ib,1)
 257            if (umbrella(iparm)) then
 258              nRR(ib,iparm)=1
 259            else
 260              nRR(ib,iparm)=nR(ib,1)
 261            endif
 262            beta_h(ib,iparm)=beta_h(ib,1)
 263            do i=1,nR(ib,iparm)
 264              f(i,ib,iparm)=f(i,ib,1)
 265              do j=1,nQ
 266                KH(j,i,ib,iparm)=KH(j,i,ib,1)
 267                Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
 268              enddo
 269            enddo
 270            replica(iparm)=replica(1)
 271            umbrella(iparm)=umbrella(1)
 272            read_iset(iparm)=read_iset(1)
 273          enddo
 274        enddo
 275
 276       endif
 277
 278       return
 279       end
 280 c-----------------------------------------------------------------------------
 281       subroutine read_protein_data(*)
 282       implicit none
 283       include "DIMENSIONS"
 284       include "DIMENSIONS.ZSCOPT"
 285       include "DIMENSIONS.FREE"
 286 #ifdef MPI
 287       include "mpif.h"
 288       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 289       include "COMMON.MPI"
 290 #endif
 291       include "COMMON.CHAIN"
 292       include "COMMON.IOUNITS"
 293       include "COMMON.PROT"
 294       include "COMMON.PROTFILES"
 295       include "COMMON.NAMES"
 296       include "COMMON.FREE"
 297       include "COMMON.OBCINKA"
 298       character*64 nazwa
 299       character*16000 controlcard
 300       integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
 301       external ilen,iroof
 302       if (hamil_rep) then
 303         npars=1
 304       else
 305         npars=nparmset
 306       endif
 307
 308       do iparm=1,npars
 309
 310 C Read names of files with conformation data.
 311       if (replica(iparm)) then
 312         nthr = 1
 313       else
 314         nthr = nT_h(iparm)
 315       endif
 316       do ib=1,nthr
 317       do ii=1,nRR(ib,iparm)
 318       write (iout,*) "Parameter set",iparm," temperature",ib,
 319      & " window",ii
 320       call flush(iout)
 321       call card_concat(controlcard,.true.)
 322       write (iout,*) controlcard(:ilen(controlcard))
 323       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
 324       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
 325       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
 326       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
 327       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
 328      & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
 329       call reada(controlcard,"TIME_START",
 330      &  time_start_collect(ii,ib,iparm),0.0d0)
 331       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
 332      &  1.0d10)
 333       write (iout,*) "rec_start",rec_start(ii,ib,iparm),
 334      & " rec_end",rec_end(ii,ib,iparm)
 335       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
 336      & " time_end",time_end_collect(ii,ib,iparm)
 337       call flush(iout)
 338       if (replica(iparm)) then
 339         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
 340         write (iout,*) "Number of trajectories",totraj(ii,iparm)
 341         call flush(iout)
 342       endif
 343       if (nfile_bin(ii,ib,iparm).lt.2
 344      &    .and. nfile_asc(ii,ib,iparm).eq.0
 345      &    .and. nfile_cx(ii,ib,iparm).eq.0) then
 346         write (iout,*) "Error - no action specified!"
 347         return1
 348       endif
 349       if (nfile_bin(ii,ib,iparm).gt.0) then
 350         call card_concat(controlcard,.false.)
 351         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
 352      &   maxfile_prot,nfile_bin(ii,ib,iparm))
 353 #ifdef DEBUG
 354         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
 355         write(iout,*) (protfiles(i,1,ii,ib,iparm),
 356      &    i=1,nfile_bin(ii,ib,iparm))
 357 #endif
 358       endif
 359       if (nfile_asc(ii,ib,iparm).gt.0) then
 360         call card_concat(controlcard,.false.)
 361         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
 362      &   maxfile_prot,nfile_asc(ii,ib,iparm))
 363 #ifdef DEBUG
 364         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
 365         write(iout,*) (protfiles(i,2,ii,ib,iparm),
 366      &    i=1,nfile_asc(ii,ib,iparm))
 367 #endif
 368       else if (nfile_cx(ii,ib,iparm).gt.0) then
 369         call card_concat(controlcard,.false.)
 370         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
 371      &   maxfile_prot,nfile_cx(ii,ib,iparm))
 372 #ifdef DEBUG
 373         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
 374         write(iout,*) (protfiles(i,2,ii,ib,iparm),
 375      &    i=1,nfile_cx(ii,ib,iparm))
 376 #endif
 377       endif
 378       call flush(iout)
 379       enddo
 380       enddo
 381
 382       enddo
 383
 384       return
 385       end
 386 c-------------------------------------------------------------------------------
 387       subroutine opentmp(islice,iunit,bprotfile_temp)
 388       implicit none
 389       include "DIMENSIONS"
 390       include "DIMENSIONS.ZSCOPT"
 391       include "DIMENSIONS.FREE"
 392 #ifdef MPI
 393       include "mpif.h"
 394       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 395       include "COMMON.MPI"
 396 #endif
 397       include "COMMON.IOUNITS"
 398       include "COMMON.PROTFILES"
 399       include "COMMON.PROT"
 400       include "COMMON.FREE"
 401       character*64 bprotfile_temp
 402       character*3 liczba,liczba2
 403       character*2 liczba1
 404       integer iunit,islice
 405       integer ilen,iroof
 406       external ilen,iroof
 407       logical lerr
 408
 409       write (liczba1,'(bz,i2.2)') islice
 410       write (liczba,'(bz,i3.3)') me
 411 #ifdef MPI
 412 c      write (iout,*) "separate_parset ",separate_parset,
 413 c     &  " myparm",myparm
 414       if (separate_parset) then
 415       write (liczba2,'(bz,i3.3)') myparm
 416       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
 417      &  prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
 418       open (iunit,file=bprotfile_temp,status="unknown",
 419      &    form="unformatted",access="direct",recl=lenrec)
 420       else
 421       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
 422      &  prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
 423       open (iunit,file=bprotfile_temp,status="unknown",
 424      &    form="unformatted",access="direct",recl=lenrec)
 425       endif
 426 #else
 427       bprotfile_temp = scratchdir(:ilen(scratchdir))//
 428      &  "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
 429       open (iunit,file=bprotfile_temp,status="unknown",
 430      &    form="unformatted",access="direct",recl=lenrec)
 431 #endif
 432 c      write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
 433 c     &  bprotfile_temp
 434 c      call flush(iout)
 435       return
 436       end
 437 c-------------------------------------------------------------------------------
 438       subroutine read_database(*)
 439       implicit none
 440       include "DIMENSIONS"
 441       include "DIMENSIONS.ZSCOPT"
 442       include "DIMENSIONS.FREE"
 443 #ifdef MPI
 444       include "mpif.h"
 445       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 446       include "COMMON.MPI"
 447 #endif
 448       include "COMMON.CHAIN"
 449       include "COMMON.IOUNITS"
 450       include "COMMON.PROTFILES"
 451       include "COMMON.NAMES"
 452       include "COMMON.VAR"
 453       include "COMMON.GEO"
 454       include "COMMON.ENEPS"
 455       include "COMMON.PROT"
 456       include "COMMON.INTERACT"
 457       include "COMMON.FREE"
 458       include "COMMON.SBRIDGE"
 459       include "COMMON.OBCINKA"
 460       real*4 csingle(3,maxres2)
 461       character*64 nazwa,bprotfile_temp
 462       character*3 liczba
 463       character*2 liczba1
 464       integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
 465      &  ll(maxslice),mm(maxslice),if
 466       integer nrec,nlines,iscor,iunit,islice
 467       double precision energ
 468       integer ilen,iroof
 469       external ilen,iroof
 470       double precision rmsdev,energia(0:max_ene),efree,eini,temp
 471       double precision prop(maxQ)
 472       integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff
 473       integer iparm,ib,iib,ir,nprop,nthr,npars
 474       double precision etot,time
 475       integer ixdrf,iret
 476       logical lerr,linit
 477
 478       lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
 479       lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
 480       lenrec=lenrec2+8
 481       write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
 482      &  " lenrec2",lenrec2
 483
 484       do i=1,nQ
 485         prop(i)=0.0d0
 486       enddo
 487       do islice=1,nslice
 488         ll(islice)=0
 489         mm(islice)=0
 490       enddo
 491       write (iout,*) "nparmset",nparmset
 492       if (hamil_rep) then
 493         npars=1
 494       else
 495         npars=nparmset
 496       endif
 497       do iparm=1,npars
 498
 499       if (replica(iparm)) then
 500         nthr = 1
 501       else
 502         nthr = nT_h(iparm)
 503       endif
 504
 505       do ib=1,nthr
 506       do iR=1,nRR(ib,iparm)
 507
 508       write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 509       do islice=1,nslice
 510         jj(islice)=0
 511         kk(islice)=0
 512       enddo
 513
 514       IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
 515 c Read conformations from binary DA files (one per batch) and write them to
 516 c a binary DA scratchfile.
 517         write (liczba,'(bz,i3.3)') me
 518         do if=1,nfile_bin(iR,ib,iparm)
 519           nazwa=protfiles(if,1,iR,ib,iparm)
 520      &     (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
 521           open (ientin,file=nazwa,status="old",form="unformatted",
 522      &     access="direct",recl=lenrec2,err=1111)
 523           ii=0
 524           do islice=1,nslice
 525             call opentmp(islice,ientout,bprotfile_temp)
 526             call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
 527      &        mm(islice),iR,ib,iparm)
 528             close(ientout)
 529           enddo
 530           close(ientin)
 531         enddo
 532       ENDIF ! NFILE_BIN>0
 533 c
 534       IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
 535 c Read conformations from multiple ASCII int files and write them to a binary
 536 c DA scratchfile.
 537         do if=1,nfile_asc(iR,ib,iparm)
 538           nazwa=protfiles(if,2,iR,ib,iparm)
 539      &     (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
 540           open(unit=ientin,file=nazwa,status='old',err=1111)
 541           write(iout,*) "reading ",nazwa(:ilen(nazwa))
 542           ii=0
 543           call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
 544         enddo ! if
 545       ENDIF
 546       IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
 547 c Read conformations from cx files and write them to a binary
 548 c DA scratchfile.
 549         do if=1,nfile_cx(iR,ib,iparm)
 550           nazwa=protfiles(if,2,iR,ib,iparm)
 551      &     (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
 552           write(iout,*) "reading ",nazwa(:ilen(nazwa))
 553           ii=0
 554           print *,"Calling cxread"
 555           call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
 556      &       *1111)
 557           close(ientout)
 558           write (iout,*) "exit cxread"
 559           call flush(iout)
 560         enddo
 561       ENDIF
 562
 563       do islice=1,nslice
 564         stot(islice)=stot(islice)+jj(islice)
 565       enddo
 566
 567       enddo
 568       enddo
 569       write (iout,*) "IPARM",iparm
 570       enddo
 571
 572       if (nslice.eq.1) then
 573 #ifdef MPI
 574         write (liczba,'(bz,i3.3)') me
 575         bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
 576      &    prefix(:ilen(prefix))//liczba//".xbin.tmp"
 577 #else
 578         bprotfile_temp = scratchdir(:ilen(scratchdir))//
 579      &     "/"//prefix(:ilen(prefix))//".xbin.tmp"
 580 #endif
 581         write(iout,*) mm(1)," conformations read",ll(1),
 582      &    " conformations written to ",
 583      &    bprotfile_temp(:ilen(bprotfile_temp))
 584       else
 585         do islice=1,nslice
 586           write (liczba1,'(bz,i2.2)') islice
 587 #ifdef MPI
 588           write (liczba,'(bz,i3.3)') me
 589           bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
 590      &      prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
 591 #else
 592           bprotfile_temp = scratchdir(:ilen(scratchdir))//
 593      &       "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
 594 #endif
 595           write(iout,*) mm(islice)," conformations read",ll(islice),
 596      &    " conformations written to ",
 597      &    bprotfile_temp(:ilen(bprotfile_temp))
 598         enddo
 599       endif
 600
 601 #ifdef MPI
 602 c Check if everyone has the same number of conformations
 603
 604 c      call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL,
 605 c     &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
 606
 607       maxslice_buff=maxslice
 608
 609       call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER,
 610      &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
 611       lerr=.false.
 612       do i=0,nprocs-1
 613         if (i.ne.me) then
 614           do islice=1,nslice
 615           if (stot(islice).ne.ntot_all(islice,i)) then
 616             write (iout,*) "Number of conformations at processor",i,
 617      &       " differs from that at processor",me,
 618      &       stot(islice),ntot_all(islice,i)," slice",islice
 619             lerr = .true.
 620           endif
 621           enddo
 622         endif
 623       enddo
 624       if (lerr) then
 625         write (iout,*)
 626         write (iout,*) "Numbers of conformations read by processors"
 627         write (iout,*)
 628         do i=0,nprocs-1
 629           write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
 630         enddo
 631         write (iout,*) "Calculation terminated."
 632         call flush(iout)
 633         return1
 634       endif
 635       do islice=1,nslice
 636         ntot(islice)=stot(islice)
 637       enddo
 638       return
 639 #endif
 640  1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
 641       call flush(iout)
 642       return1
 643       end
 644 c------------------------------------------------------------------------------
 645       subroutine card_concat(card,to_upper)
 646       implicit none
 647       include 'DIMENSIONS.ZSCOPT'
 648       include "COMMON.IOUNITS"
 649       character*(*) card
 650       character*80 karta,ucase
 651       logical to_upper
 652       integer ilen
 653       external ilen
 654       read (inp,'(a)') karta
 655       if (to_upper) karta=ucase(karta)
 656       card=' '
 657       do while (karta(80:80).eq.'&')
 658         card=card(:ilen(card)+1)//karta(:79)
 659         read (inp,'(a)') karta
 660         if (to_upper) karta=ucase(karta)
 661       enddo
 662       card=card(:ilen(card)+1)//karta
 663       return
 664       end
 665 c------------------------------------------------------------------------------
 666       subroutine readi(rekord,lancuch,wartosc,default)
 667       implicit none
 668       character*(*) rekord,lancuch
 669       integer wartosc,default
 670       integer ilen,iread
 671       external ilen
 672       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 673       if (iread.eq.0) then
 674         wartosc=default
 675         return
 676       endif
 677       iread=iread+ilen(lancuch)+1
 678       read (rekord(iread:),*) wartosc
 679       return
 680       end
 681 c----------------------------------------------------------------------------
 682       subroutine reada(rekord,lancuch,wartosc,default)
 683       implicit none
 684       character*(*) rekord,lancuch
 685       character*80 aux
 686       double precision wartosc,default
 687       integer ilen,iread
 688       external ilen
 689       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 690       if (iread.eq.0) then
 691         wartosc=default
 692         return
 693       endif
 694       iread=iread+ilen(lancuch)+1
 695       read (rekord(iread:),*) wartosc
 696       return
 697       end
 698 c----------------------------------------------------------------------------
 699       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 700       implicit none
 701       integer dim,i
 702       integer tablica(dim),default
 703       character*(*) rekord,lancuch
 704       character*80 aux
 705       integer ilen,iread
 706       external ilen
 707       do i=1,dim
 708         tablica(i)=default
 709       enddo
 710       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 711       if (iread.eq.0) return
 712       iread=iread+ilen(lancuch)+1
 713       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 714    10 return
 715       end
 716 c----------------------------------------------------------------------------
 717       subroutine multreada(rekord,lancuch,tablica,dim,default)
 718       implicit none
 719       integer dim,i
 720       double precision tablica(dim),default
 721       character*(*) rekord,lancuch
 722       character*80 aux
 723       integer ilen,iread
 724       external ilen
 725       do i=1,dim
 726         tablica(i)=default
 727       enddo
 728       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 729       if (iread.eq.0) return
 730       iread=iread+ilen(lancuch)+1
 731       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 732    10 return
 733       end
 734 c----------------------------------------------------------------------------
 735       subroutine reads(rekord,lancuch,wartosc,default)
 736       implicit none
 737       character*(*) rekord,lancuch,wartosc,default
 738       character*80 aux
 739       integer ilen,lenlan,lenrec,iread,ireade
 740       external ilen
 741       logical iblnk
 742       external iblnk
 743       lenlan=ilen(lancuch)
 744       lenrec=ilen(rekord)
 745       iread=index(rekord,lancuch(:lenlan)//"=")
 746 c      print *,"rekord",rekord," lancuch",lancuch
 747 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
 748       if (iread.eq.0) then
 749         wartosc=default
 750         return
 751       endif
 752       iread=iread+lenlan+1
 753 c      print *,"iread",iread
 754 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
 755       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
 756         iread=iread+1
 757 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
 758       enddo
 759 c      print *,"iread",iread
 760       if (iread.gt.lenrec) then
 761          wartosc=default
 762         return
 763       endif
 764       ireade=iread+1
 765 c      print *,"ireade",ireade
 766       do while (ireade.lt.lenrec .and.
 767      &   .not.iblnk(rekord(ireade:ireade)))
 768         ireade=ireade+1
 769       enddo
 770       wartosc=rekord(iread:ireade)
 771       return
 772       end
 773 c----------------------------------------------------------------------------
 774       subroutine multreads(rekord,lancuch,tablica,dim,default)
 775       implicit none
 776       integer dim,i
 777       character*(*) rekord,lancuch,tablica(dim),default
 778       character*80 aux
 779       integer ilen,lenlan,lenrec,iread,ireade
 780       external ilen
 781       logical iblnk
 782       external iblnk
 783       do i=1,dim
 784         tablica(i)=default
 785       enddo
 786       lenlan=ilen(lancuch)
 787       lenrec=ilen(rekord)
 788       iread=index(rekord,lancuch(:lenlan)//"=")
 789 c      print *,"rekord",rekord," lancuch",lancuch
 790 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
 791       if (iread.eq.0) return
 792       iread=iread+lenlan+1
 793       do i=1,dim
 794 c      print *,"iread",iread
 795 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
 796       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
 797         iread=iread+1
 798 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
 799       enddo
 800 c      print *,"iread",iread
 801       if (iread.gt.lenrec) return
 802       ireade=iread+1
 803 c      print *,"ireade",ireade
 804       do while (ireade.lt.lenrec .and.
 805      &   .not.iblnk(rekord(ireade:ireade)))
 806         ireade=ireade+1
 807       enddo
 808       tablica(i)=rekord(iread:ireade)
 809       iread=ireade+1
 810       enddo
 811       end
 812 c----------------------------------------------------------------------------
 813       subroutine split_string(rekord,tablica,dim,nsub)
 814       implicit none
 815       integer dim,nsub,i,ii,ll,kk
 816       character*(*) tablica(dim)
 817       character*(*) rekord
 818       integer ilen
 819       external ilen
 820       do i=1,dim
 821         tablica(i)=" "
 822       enddo
 823       ii=1
 824       ll = ilen(rekord)
 825       nsub=0
 826       do i=1,dim
 827 C Find the start of term name
 828         kk = 0
 829         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
 830           ii = ii+1
 831         enddo
 832 C Parse the name into TABLICA(i) until blank found
 833         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
 834           kk = kk+1
 835           tablica(i)(kk:kk)=rekord(ii:ii)
 836           ii = ii+1
 837         enddo
 838         if (kk.gt.0) nsub=nsub+1
 839         if (ii.gt.ll) return
 840       enddo
 841       return
 842       end
 843 c--------------------------------------------------------------------------------
 844       integer function iroof(n,m)
 845       ii = n/m
 846       if (ii*m .lt. n) ii=ii+1
 847       iroof = ii
 848       return
 849       end