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