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