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