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