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