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