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