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