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