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       do i=1,nparmset
  55         if (isampl(i).eq.0) then
  56           write (iout,*) "ERROR: isampl is 0 for parmset",i
  57           call flush(iout)
  58           stop
  59         endif
  60       enddo
  61       write (iout,*) "MaxSlice",MaxSlice
  62       call readi(controlcard,"NSLICE",nslice,1)
  63       call flush(iout)
  64       if (nslice.gt.MaxSlice) then
  65         write (iout,*) "Error: parameter out of range: NSLICE",nslice,
  66      &    MaxSlice
  67         return1
  68       endif
  69       write (iout,*) "Frequency of storing conformations",
  70      & (isampl(i),i=1,nparmset)
  71       write (iout,*) "Maxit",maxit," Fimin",fimin
  72       call readi(controlcard,"NQ",nQ,1)
  73       if (nQ.gt.MaxQ) then
  74         write (iout,*) "Error: parameter out of range: NQ",nq,
  75      &    maxq
  76         return1
  77       endif
  78       indpdb=0
  79       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
  80       rmsrgymap = (index(controlcard,'RMSRGYMAP').gt.0)
  81       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
  82       call reada(controlcard,"DELTA",delta,1.0d-2)
  83       call readi(controlcard,"EINICHECK",einicheck,2)
  84       call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
  85       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
  86       call readi(controlcard,"RESCALE",rescale_mode,1)
  87       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
  88       call readi(controlcard,'TORMODE',tor_mode,0)
  89 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
  90       write(iout,*) "torsional and valence angle mode",tor_mode
  91       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
  92        call reada(controlcard,'BOXX',boxxsize,100.0d0)
  93        call reada(controlcard,'BOXY',boxysize,100.0d0)
  94        call reada(controlcard,'BOXZ',boxzsize,100.0d0)
  95 c Cutoff range for interactions
  96       call reada(controlcard,"R_CUT",r_cut,15.0d0)
  97       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
  98       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
  99       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 100       if (lipthick.gt.0.0d0) then
 101        bordliptop=(boxzsize+lipthick)/2.0
 102        bordlipbot=bordliptop-lipthick
 103 C      endif
 104       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 105      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 106       buflipbot=bordlipbot+lipbufthick
 107       bufliptop=bordliptop-lipbufthick
 108       if ((lipbufthick*2.0d0).gt.lipthick)
 109      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 110       endif
 111       write(iout,*) "bordliptop=",bordliptop
 112       write(iout,*) "bordlipbot=",bordlipbot
 113       write(iout,*) "bufliptop=",bufliptop
 114       write(iout,*) "buflipbot=",buflipbot
 115       call readi(controlcard,'SYM',symetr,1)
 116       write (iout,*) "DISTCHAINMAX",distchainmax
 117       write (iout,*) "delta",delta
 118       write (iout,*) "einicheck",einicheck
 119       write (iout,*) "rescale_mode",rescale_mode
 120       call flush(iout)
 121       bxfile=index(controlcard,"BXFILE").gt.0
 122       cxfile=index(controlcard,"CXFILE").gt.0
 123       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
 124      & bxfile=.true.
 125       histfile=index(controlcard,"HISTFILE").gt.0
 126       histout=index(controlcard,"HISTOUT").gt.0
 127       entfile=index(controlcard,"ENTFILE").gt.0
 128       zscfile=index(controlcard,"ZSCFILE").gt.0
 129       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
 130       write (iout,*) "with_dihed_constr ",with_dihed_constr
 131       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 132       write (iout,*) "with_theta_constr ",with_theta_constr
 133       call readi(controlcard,'SHIELD',shield_mode,0)
 134         write(iout,*) "shield_mode",shield_mode
 135 C      endif
 136       call readi(controlcard,'TORMODE',tor_mode,0)
 137         write(iout,*) "torsional and valence angle mode",tor_mode
 138       if (shield_mode.gt.0) then
 139       pi=3.141592d0
 140 C VSolvSphere the volume of solving sphere
 141 C      print *,pi,"pi"
 142 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 143 C there will be no distinction between proline peptide group and normal peptide
 144 C group in case of shielding parameters
 145       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 146       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 147       write (iout,*) VSolvSphere,VSolvSphere_div
 148 C long axis of side chain
 149 C      do i=1,ntyp
 150 C      long_r_sidechain(i)=vbldsc0(1,i)
 151 C      short_r_sidechain(i)=sigma0(i)
 152 C      enddo
 153       buff_shield=1.0d0
 154       endif
 155
 156       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 157       dyn_ss=index(controlcard,"DYN_SS").gt.0
 158       adaptive = index(controlcard,"ADAPTIVE").gt.0
 159       return
 160       end
 161 c------------------------------------------------------------------------------
 162       subroutine read_efree(*)
 163 C
 164 C Read molecular data
 165 C
 166       implicit none
 167       include 'DIMENSIONS'
 168       include 'DIMENSIONS.ZSCOPT'
 169       include 'DIMENSIONS.COMPAR'
 170       include 'DIMENSIONS.FREE'
 171       include 'COMMON.IOUNITS'
 172       include 'COMMON.TIME1'
 173       include 'COMMON.SBRIDGE'
 174       include 'COMMON.CONTROL'
 175       include 'COMMON.CHAIN'
 176       include 'COMMON.HEADER'
 177       include 'COMMON.GEO'
 178       include 'COMMON.FREE'
 179       character*320 controlcard,ucase
 180       integer iparm,ib,i,j,npars
 181       integer ilen
 182       external ilen
 183
 184       if (hamil_rep) then
 185         npars=1
 186       else
 187         npars=nParmSet
 188       endif
 189
 190       do iparm=1,npars
 191
 192       call card_concat(controlcard,.true.)
 193       call readi(controlcard,'NT',nT_h(iparm),1)
 194       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
 195       call flush(iout)
 196       if (nT_h(iparm).gt.MaxT_h) then
 197         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),
 198      &    MaxT_h
 199         return1
 200       endif
 201       replica(iparm)=index(controlcard,"REPLICA").gt.0
 202       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
 203       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
 204       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
 205      &  replica(iparm)," umbrella ",umbrella(iparm),
 206      &  " read_iset",read_iset(iparm)
 207       call flush(iout)
 208       do ib=1,nT_h(iparm)
 209         call card_concat(controlcard,.true.)
 210         call readi(controlcard,'NR',nR(ib,iparm),1)
 211         if (umbrella(iparm)) then
 212           nRR(ib,iparm)=1
 213         else
 214           nRR(ib,iparm)=nR(ib,iparm)
 215         endif
 216         if (nR(ib,iparm).gt.MaxR) then
 217           write (iout,*)  "Error: parameter out of range: NR",
 218      &      nR(ib,iparm),MaxR
 219         return1
 220         endif
 221         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
 222         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
 223         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
 224      &    0.0d0)
 225         do i=1,nR(ib,iparm)
 226           call card_concat(controlcard,.true.)
 227           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
 228      &      100.0d0)
 229           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
 230      &      0.0d0)
 231         enddo
 232       enddo
 233       do ib=1,nT_h(iparm)
 234         write (iout,*) "ib",ib," beta_h",
 235      &    1.0d0/(0.001987*beta_h(ib,iparm))
 236         write (iout,*) "nR",nR(ib,iparm)
 237         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
 238         do i=1,nR(ib,iparm)
 239           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
 240      &      "q0",(q0(j,i,ib,iparm),j=1,nQ)
 241         enddo
 242         call flush(iout)
 243       enddo
 244
 245       enddo
 246
 247       if (hamil_rep) then
 248
 249        do iparm=2,nParmSet
 250           nT_h(iparm)=nT_h(1)
 251          do ib=1,nT_h(iparm)
 252            nR(ib,iparm)=nR(ib,1)
 253            if (umbrella(iparm)) then
 254              nRR(ib,iparm)=1
 255            else
 256              nRR(ib,iparm)=nR(ib,1)
 257            endif
 258            beta_h(ib,iparm)=beta_h(ib,1)
 259            do i=1,nR(ib,iparm)
 260              f(i,ib,iparm)=f(i,ib,1)
 261              do j=1,nQ
 262                KH(j,i,ib,iparm)=KH(j,i,ib,1)
 263                Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
 264              enddo
 265            enddo
 266            replica(iparm)=replica(1)
 267            umbrella(iparm)=umbrella(1)
 268            read_iset(iparm)=read_iset(1)
 269          enddo
 270        enddo
 271
 272       endif
 273
 274       return
 275       end
 276 c-----------------------------------------------------------------------------
 277       subroutine read_protein_data(*)
 278       implicit none
 279       include "DIMENSIONS"
 280       include "DIMENSIONS.ZSCOPT"
 281       include "DIMENSIONS.FREE"
 282 #ifdef MPI
 283       include "mpif.h"
 284       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 285       include "COMMON.MPI"
 286 #endif
 287       include "COMMON.CHAIN"
 288       include "COMMON.IOUNITS"
 289       include "COMMON.PROT"
 290       include "COMMON.PROTFILES"
 291       include "COMMON.NAMES"
 292       include "COMMON.FREE"
 293       include "COMMON.OBCINKA"
 294       character*64 nazwa
 295       character*16000 controlcard
 296       integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
 297       external ilen,iroof
 298       if (hamil_rep) then
 299         npars=1
 300       else
 301         npars=nparmset
 302       endif
 303
 304       do iparm=1,npars
 305
 306 C Read names of files with conformation data.
 307       if (replica(iparm)) then
 308         nthr = 1
 309       else
 310         nthr = nT_h(iparm)
 311       endif
 312       do ib=1,nthr
 313       do ii=1,nRR(ib,iparm)
 314       write (iout,*) "Parameter set",iparm," temperature",ib,
 315      & " window",ii
 316       call flush(iout)
 317       call card_concat(controlcard,.true.)
 318       write (iout,*) controlcard(:ilen(controlcard))
 319       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
 320       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
 321       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
 322       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
 323       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
 324      & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
 325       call reada(controlcard,"TIME_START",
 326      &  time_start_collect(ii,ib,iparm),0.0d0)
 327       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
 328      &  1.0d10)
 329       write (iout,*) "rec_start",rec_start(ii,ib,iparm),
 330      & " rec_end",rec_end(ii,ib,iparm)
 331       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
 332      & " time_end",time_end_collect(ii,ib,iparm)
 333       call flush(iout)
 334       if (replica(iparm)) then
 335         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
 336         write (iout,*) "Number of trajectories",totraj(ii,iparm)
 337         call flush(iout)
 338       endif
 339       if (nfile_bin(ii,ib,iparm).lt.2
 340      &    .and. nfile_asc(ii,ib,iparm).eq.0
 341      &    .and. nfile_cx(ii,ib,iparm).eq.0) then
 342         write (iout,*) "Error - no action specified!"
 343         return1
 344       endif
 345       if (nfile_bin(ii,ib,iparm).gt.0) then
 346         call card_concat(controlcard,.false.)
 347         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
 348      &   maxfile_prot,nfile_bin(ii,ib,iparm))
 349 #ifdef DEBUG
 350         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
 351         write(iout,*) (protfiles(i,1,ii,ib,iparm),
 352      &    i=1,nfile_bin(ii,ib,iparm))
 353 #endif
 354       endif
 355       if (nfile_asc(ii,ib,iparm).gt.0) then
 356         call card_concat(controlcard,.false.)
 357         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
 358      &   maxfile_prot,nfile_asc(ii,ib,iparm))
 359 #ifdef DEBUG
 360         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
 361         write(iout,*) (protfiles(i,2,ii,ib,iparm),
 362      &    i=1,nfile_asc(ii,ib,iparm))
 363 #endif
 364       else if (nfile_cx(ii,ib,iparm).gt.0) then
 365         call card_concat(controlcard,.false.)
 366         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
 367      &   maxfile_prot,nfile_cx(ii,ib,iparm))
 368 #ifdef DEBUG
 369         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
 370         write(iout,*) (protfiles(i,2,ii,ib,iparm),
 371      &    i=1,nfile_cx(ii,ib,iparm))
 372 #endif
 373       endif
 374       call flush(iout)
 375       enddo
 376       enddo
 377
 378       enddo
 379
 380       return
 381       end
 382 c-------------------------------------------------------------------------------
 383       subroutine opentmp(islice,iunit,bprotfile_temp)
 384       implicit none
 385       include "DIMENSIONS"
 386       include "DIMENSIONS.ZSCOPT"
 387       include "DIMENSIONS.FREE"
 388 #ifdef MPI
 389       include "mpif.h"
 390       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 391       include "COMMON.MPI"
 392 #endif
 393       include "COMMON.IOUNITS"
 394       include "COMMON.PROTFILES"
 395       include "COMMON.PROT"
 396       include "COMMON.FREE"
 397       character*64 bprotfile_temp
 398       character*3 liczba,liczba2
 399       character*2 liczba1
 400       integer iunit,islice
 401       integer ilen,iroof
 402       external ilen,iroof
 403       logical lerr
 404
 405       write (liczba1,'(bz,i2.2)') islice
 406       write (liczba,'(bz,i3.3)') me
 407 #ifdef MPI
 408 c      write (iout,*) "separate_parset ",separate_parset,
 409 c     &  " myparm",myparm
 410       if (separate_parset) then
 411       write (liczba2,'(bz,i3.3)') myparm
 412       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
 413      &  prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
 414       open (iunit,file=bprotfile_temp,status="unknown",
 415      &    form="unformatted",access="direct",recl=lenrec)
 416       else
 417       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
 418      &  prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
 419       open (iunit,file=bprotfile_temp,status="unknown",
 420      &    form="unformatted",access="direct",recl=lenrec)
 421       endif
 422 #else
 423       bprotfile_temp = scratchdir(:ilen(scratchdir))//
 424      &  "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
 425       open (iunit,file=bprotfile_temp,status="unknown",
 426      &    form="unformatted",access="direct",recl=lenrec)
 427 #endif
 428 c      write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
 429 c     &  bprotfile_temp
 430 c      call flush(iout)
 431       return
 432       end
 433 c-------------------------------------------------------------------------------
 434       subroutine read_database(*)
 435       implicit none
 436       include "DIMENSIONS"
 437       include "DIMENSIONS.ZSCOPT"
 438       include "DIMENSIONS.FREE"
 439 #ifdef MPI
 440       include "mpif.h"
 441       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 442       include "COMMON.MPI"
 443 #endif
 444       include "COMMON.CHAIN"
 445       include "COMMON.IOUNITS"
 446       include "COMMON.PROTFILES"
 447       include "COMMON.NAMES"
 448       include "COMMON.VAR"
 449       include "COMMON.GEO"
 450       include "COMMON.ENEPS"
 451       include "COMMON.PROT"
 452       include "COMMON.INTERACT"
 453       include "COMMON.FREE"
 454       include "COMMON.SBRIDGE"
 455       include "COMMON.OBCINKA"
 456       real*4 csingle(3,maxres2)
 457       character*64 nazwa,bprotfile_temp
 458       character*3 liczba
 459       character*2 liczba1
 460       integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
 461      &  ll(maxslice),mm(maxslice),if
 462       integer nrec,nlines,iscor,iunit,islice
 463       double precision energ
 464       integer ilen,iroof
 465       external ilen,iroof
 466       double precision rmsdev,energia(0:max_ene),efree,eini,temp
 467       double precision prop(maxQ)
 468       integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff
 469       integer iparm,ib,iib,ir,nprop,nthr,npars
 470       double precision etot,time
 471       integer ixdrf,iret
 472       logical lerr,linit
 473
 474       lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
 475       lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
 476       lenrec=lenrec2+8
 477       write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
 478      &  " lenrec2",lenrec2
 479
 480       do i=1,nQ
 481         prop(i)=0.0d0
 482       enddo
 483       do islice=1,nslice
 484         ll(islice)=0
 485         mm(islice)=0
 486       enddo
 487       write (iout,*) "nparmset",nparmset
 488       if (hamil_rep) then
 489         npars=1
 490       else
 491         npars=nparmset
 492       endif
 493       do iparm=1,npars
 494
 495       if (replica(iparm)) then
 496         nthr = 1
 497       else
 498         nthr = nT_h(iparm)
 499       endif
 500
 501       do ib=1,nthr
 502       do iR=1,nRR(ib,iparm)
 503
 504       write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 505       do islice=1,nslice
 506         jj(islice)=0
 507         kk(islice)=0
 508       enddo
 509
 510       IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
 511 c Read conformations from binary DA files (one per batch) and write them to
 512 c a binary DA scratchfile.
 513         write (liczba,'(bz,i3.3)') me
 514         do if=1,nfile_bin(iR,ib,iparm)
 515           nazwa=protfiles(if,1,iR,ib,iparm)
 516      &     (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
 517           open (ientin,file=nazwa,status="old",form="unformatted",
 518      &     access="direct",recl=lenrec2,err=1111)
 519           ii=0
 520           do islice=1,nslice
 521             call opentmp(islice,ientout,bprotfile_temp)
 522             call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
 523      &        mm(islice),iR,ib,iparm)
 524             close(ientout)
 525           enddo
 526           close(ientin)
 527         enddo
 528       ENDIF ! NFILE_BIN>0
 529 c
 530       IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
 531 c Read conformations from multiple ASCII int files and write them to a binary
 532 c DA scratchfile.
 533         do if=1,nfile_asc(iR,ib,iparm)
 534           nazwa=protfiles(if,2,iR,ib,iparm)
 535      &     (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
 536           open(unit=ientin,file=nazwa,status='old',err=1111)
 537           write(iout,*) "reading ",nazwa(:ilen(nazwa))
 538           ii=0
 539           call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
 540         enddo ! if
 541       ENDIF
 542       IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
 543 c Read conformations from cx files and write them to a binary
 544 c DA scratchfile.
 545         do if=1,nfile_cx(iR,ib,iparm)
 546           nazwa=protfiles(if,2,iR,ib,iparm)
 547      &     (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
 548           write(iout,*) "reading ",nazwa(:ilen(nazwa))
 549           ii=0
 550           print *,"Calling cxread"
 551           call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
 552      &       *1111)
 553           close(ientout)
 554           write (iout,*) "exit cxread"
 555           call flush(iout)
 556         enddo
 557       ENDIF
 558
 559       do islice=1,nslice
 560         stot(islice)=stot(islice)+jj(islice)
 561       enddo
 562
 563       enddo
 564       enddo
 565       write (iout,*) "IPARM",iparm
 566       enddo
 567
 568       if (nslice.eq.1) then
 569 #ifdef MPI
 570         write (liczba,'(bz,i3.3)') me
 571         bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
 572      &    prefix(:ilen(prefix))//liczba//".xbin.tmp"
 573 #else
 574         bprotfile_temp = scratchdir(:ilen(scratchdir))//
 575      &     "/"//prefix(:ilen(prefix))//".xbin.tmp"
 576 #endif
 577         write(iout,*) mm(1)," conformations read",ll(1),
 578      &    " conformations written to ",
 579      &    bprotfile_temp(:ilen(bprotfile_temp))
 580       else
 581         do islice=1,nslice
 582           write (liczba1,'(bz,i2.2)') islice
 583 #ifdef MPI
 584           write (liczba,'(bz,i3.3)') me
 585           bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
 586      &      prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
 587 #else
 588           bprotfile_temp = scratchdir(:ilen(scratchdir))//
 589      &       "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
 590 #endif
 591           write(iout,*) mm(islice)," conformations read",ll(islice),
 592      &    " conformations written to ",
 593      &    bprotfile_temp(:ilen(bprotfile_temp))
 594         enddo
 595       endif
 596
 597 #ifdef MPI
 598 c Check if everyone has the same number of conformations
 599
 600 c      call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL,
 601 c     &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
 602
 603       maxslice_buff=maxslice
 604
 605       call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER,
 606      &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
 607       lerr=.false.
 608       do i=0,nprocs-1
 609         if (i.ne.me) then
 610           do islice=1,nslice
 611           if (stot(islice).ne.ntot_all(islice,i)) then
 612             write (iout,*) "Number of conformations at processor",i,
 613      &       " differs from that at processor",me,
 614      &       stot(islice),ntot_all(islice,i)," slice",islice
 615             lerr = .true.
 616           endif
 617           enddo
 618         endif
 619       enddo
 620       if (lerr) then
 621         write (iout,*)
 622         write (iout,*) "Numbers of conformations read by processors"
 623         write (iout,*)
 624         do i=0,nprocs-1
 625           write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
 626         enddo
 627         write (iout,*) "Calculation terminated."
 628         call flush(iout)
 629         return1
 630       endif
 631       do islice=1,nslice
 632         ntot(islice)=stot(islice)
 633       enddo
 634       return
 635 #endif
 636  1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
 637       call flush(iout)
 638       return1
 639       end
 640 c------------------------------------------------------------------------------
 641       subroutine card_concat(card,to_upper)
 642       implicit none
 643       include 'DIMENSIONS.ZSCOPT'
 644       include "COMMON.IOUNITS"
 645       character*(*) card
 646       character*80 karta,ucase
 647       logical to_upper
 648       integer ilen
 649       external ilen
 650       read (inp,'(a)') karta
 651       if (to_upper) karta=ucase(karta)
 652       card=' '
 653       do while (karta(80:80).eq.'&')
 654         card=card(:ilen(card)+1)//karta(:79)
 655         read (inp,'(a)') karta
 656         if (to_upper) karta=ucase(karta)
 657       enddo
 658       card=card(:ilen(card)+1)//karta
 659       return
 660       end
 661 c------------------------------------------------------------------------------
 662       subroutine readi(rekord,lancuch,wartosc,default)
 663       implicit none
 664       character*(*) rekord,lancuch
 665       integer 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 reada(rekord,lancuch,wartosc,default)
 679       implicit none
 680       character*(*) rekord,lancuch
 681       character*80 aux
 682       double precision wartosc,default
 683       integer ilen,iread
 684       external ilen
 685       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 686       if (iread.eq.0) then
 687         wartosc=default
 688         return
 689       endif
 690       iread=iread+ilen(lancuch)+1
 691       read (rekord(iread:),*) wartosc
 692       return
 693       end
 694 c----------------------------------------------------------------------------
 695       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 696       implicit none
 697       integer dim,i
 698       integer tablica(dim),default
 699       character*(*) rekord,lancuch
 700       character*80 aux
 701       integer ilen,iread
 702       external ilen
 703       do i=1,dim
 704         tablica(i)=default
 705       enddo
 706       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 707       if (iread.eq.0) return
 708       iread=iread+ilen(lancuch)+1
 709       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 710    10 return
 711       end
 712 c----------------------------------------------------------------------------
 713       subroutine multreada(rekord,lancuch,tablica,dim,default)
 714       implicit none
 715       integer dim,i
 716       double precision tablica(dim),default
 717       character*(*) rekord,lancuch
 718       character*80 aux
 719       integer ilen,iread
 720       external ilen
 721       do i=1,dim
 722         tablica(i)=default
 723       enddo
 724       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 725       if (iread.eq.0) return
 726       iread=iread+ilen(lancuch)+1
 727       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 728    10 return
 729       end
 730 c----------------------------------------------------------------------------
 731       subroutine reads(rekord,lancuch,wartosc,default)
 732       implicit none
 733       character*(*) rekord,lancuch,wartosc,default
 734       character*80 aux
 735       integer ilen,lenlan,lenrec,iread,ireade
 736       external ilen
 737       logical iblnk
 738       external iblnk
 739       lenlan=ilen(lancuch)
 740       lenrec=ilen(rekord)
 741       iread=index(rekord,lancuch(:lenlan)//"=")
 742 c      print *,"rekord",rekord," lancuch",lancuch
 743 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
 744       if (iread.eq.0) then
 745         wartosc=default
 746         return
 747       endif
 748       iread=iread+lenlan+1
 749 c      print *,"iread",iread
 750 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
 751       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
 752         iread=iread+1
 753 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
 754       enddo
 755 c      print *,"iread",iread
 756       if (iread.gt.lenrec) then
 757          wartosc=default
 758         return
 759       endif
 760       ireade=iread+1
 761 c      print *,"ireade",ireade
 762       do while (ireade.lt.lenrec .and.
 763      &   .not.iblnk(rekord(ireade:ireade)))
 764         ireade=ireade+1
 765       enddo
 766       wartosc=rekord(iread:ireade)
 767       return
 768       end
 769 c----------------------------------------------------------------------------
 770       subroutine multreads(rekord,lancuch,tablica,dim,default)
 771       implicit none
 772       integer dim,i
 773       character*(*) rekord,lancuch,tablica(dim),default
 774       character*80 aux
 775       integer ilen,lenlan,lenrec,iread,ireade
 776       external ilen
 777       logical iblnk
 778       external iblnk
 779       do i=1,dim
 780         tablica(i)=default
 781       enddo
 782       lenlan=ilen(lancuch)
 783       lenrec=ilen(rekord)
 784       iread=index(rekord,lancuch(:lenlan)//"=")
 785 c      print *,"rekord",rekord," lancuch",lancuch
 786 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
 787       if (iread.eq.0) return
 788       iread=iread+lenlan+1
 789       do i=1,dim
 790 c      print *,"iread",iread
 791 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
 792       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
 793         iread=iread+1
 794 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
 795       enddo
 796 c      print *,"iread",iread
 797       if (iread.gt.lenrec) return
 798       ireade=iread+1
 799 c      print *,"ireade",ireade
 800       do while (ireade.lt.lenrec .and.
 801      &   .not.iblnk(rekord(ireade:ireade)))
 802         ireade=ireade+1
 803       enddo
 804       tablica(i)=rekord(iread:ireade)
 805       iread=ireade+1
 806       enddo
 807       end
 808 c----------------------------------------------------------------------------
 809       subroutine split_string(rekord,tablica,dim,nsub)
 810       implicit none
 811       integer dim,nsub,i,ii,ll,kk
 812       character*(*) tablica(dim)
 813       character*(*) rekord
 814       integer ilen
 815       external ilen
 816       do i=1,dim
 817         tablica(i)=" "
 818       enddo
 819       ii=1
 820       ll = ilen(rekord)
 821       nsub=0
 822       do i=1,dim
 823 C Find the start of term name
 824         kk = 0
 825         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
 826           ii = ii+1
 827         enddo
 828 C Parse the name into TABLICA(i) until blank found
 829         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
 830           kk = kk+1
 831           tablica(i)(kk:kk)=rekord(ii:ii)
 832           ii = ii+1
 833         enddo
 834         if (kk.gt.0) nsub=nsub+1
 835         if (ii.gt.ll) return
 836       enddo
 837       return
 838       end
 839 c--------------------------------------------------------------------------------
 840       integer function iroof(n,m)
 841       ii = n/m
 842       if (ii*m .lt. n) ii=ii+1
 843       iroof = ii
 844       return
 845       end
 846 c-------------------------------------------------------------------------------
 847       subroutine read_dist_constr
 848       implicit real*8 (a-h,o-z)
 849       include 'DIMENSIONS'
 850       include 'COMMON.CONTROL'
 851       include 'COMMON.CHAIN'
 852       include 'COMMON.IOUNITS'
 853       include 'COMMON.SBRIDGE'
 854       integer ifrag_(2,100),ipair_(2,100)
 855       double precision wfrag_(100),wpair_(100)
 856       character*500 controlcard
 857       logical normalize
 858       print *, "WCHODZE"
 859       write (iout,*) "Calling read_dist_constr"
 860 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 861 c      call flush(iout)
 862       call card_concat(controlcard,.true.)
 863       call readi(controlcard,"NFRAG",nfrag_,0)
 864       call readi(controlcard,"NPAIR",npair_,0)
 865       call readi(controlcard,"NDIST",ndist_,0)
 866       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 867       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 868       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 869       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 870       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 871       normalize = index(controlcard,"NORMALIZE").gt.0
 872 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 873 c      write (iout,*) "IFRAG"
 874 c      do i=1,nfrag_
 875 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 876 c      enddo
 877 c      write (iout,*) "IPAIR"
 878 c      do i=1,npair_
 879 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 880 c      enddo
 881       call flush(iout)
 882       do i=1,nfrag_
 883         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 884         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 885      &    ifrag_(2,i)=nstart_sup+nsup-1
 886 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 887         call flush(iout)
 888         if (wfrag_(i).gt.0.0d0) then
 889         do j=ifrag_(1,i),ifrag_(2,i)-1
 890           do k=j+1,ifrag_(2,i)
 891 c            write (iout,*) "j",j," k",k
 892             ddjk=dist(j,k)
 893             if (constr_dist.eq.1) then
 894             nhpb=nhpb+1
 895             ihpb(nhpb)=j
 896             jhpb(nhpb)=k
 897               dhpb(nhpb)=ddjk
 898             forcon(nhpb)=wfrag_(i)
 899             else if (constr_dist.eq.2) then
 900               if (ddjk.le.dist_cut) then
 901                 nhpb=nhpb+1
 902                 ihpb(nhpb)=j
 903                 jhpb(nhpb)=k
 904                 dhpb(nhpb)=ddjk
 905                 forcon(nhpb)=wfrag_(i)
 906               endif
 907             else
 908               nhpb=nhpb+1
 909               ihpb(nhpb)=j
 910               jhpb(nhpb)=k
 911               dhpb(nhpb)=ddjk
 912               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 913             endif
 914             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 915      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 916           enddo
 917         enddo
 918         endif
 919       enddo
 920       do i=1,npair_
 921         if (wpair_(i).gt.0.0d0) then
 922         ii = ipair_(1,i)
 923         jj = ipair_(2,i)
 924         if (ii.gt.jj) then
 925           itemp=ii
 926           ii=jj
 927           jj=itemp
 928         endif
 929         do j=ifrag_(1,ii),ifrag_(2,ii)
 930           do k=ifrag_(1,jj),ifrag_(2,jj)
 931             nhpb=nhpb+1
 932             ihpb(nhpb)=j
 933             jhpb(nhpb)=k
 934             forcon(nhpb)=wpair_(i)
 935             dhpb(nhpb)=dist(j,k)
 936             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 937      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 938           enddo
 939         enddo
 940         endif
 941       enddo
 942       write (iout,*) "Distance restraints as read from input"
 943       do i=1,ndist_
 944         if (constr_dist.eq.11) then
 945           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 946      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 947 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 948           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
 949           nhpb=nhpb+1
 950           write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ",
 951      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 952      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb)
 953           if (ibecarb(i).gt.0) then
 954             ihpb(i)=ihpb(i)+nres
 955             jhpb(i)=jhpb(i)+nres
 956           endif
 957         else
 958 C        print *,"in else"
 959           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 960      &     ibecarb(i),forcon(nhpb+1)
 961           if (forcon(nhpb+1).gt.0.0d0) then
 962           nhpb=nhpb+1
 963           if (ibecarb(i).gt.0) then
 964             ihpb(i)=ihpb(i)+nres
 965             jhpb(i)=jhpb(i)+nres
 966           endif
 967           if (dhpb(nhpb).eq.0.0d0)
 968      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 969         endif
 970           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
 971      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
 972         endif
 973 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 974 C        if (forcon(nhpb+1).gt.0.0d0) then
 975 C          nhpb=nhpb+1
 976 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 977       enddo
 978       if (constr_dist.eq.11 .and. normalize) then
 979         fordepthmax=fordepth(1)
 980         do i=2,nhpb
 981           if (fordepth(i).gt.fordepthmax) fordepthmax=fordepth(i)
 982         enddo
 983         do i=1,nhpb
 984           fordepth(i)=fordepth(i)/fordepthmax
 985         enddo
 986         write (iout,'(/a/4a5,2a8,2a10)')
 987      &   "Normalized Lorenzian-like distance restraints",
 988      &   "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V"
 989         do i=1,nhpb
 990           write (iout,'(4i5,2f8.2,2f10.5)')i,ihpb(i),jhpb(i),ibecarb(i),
 991      &      dhpb(i),dhpb1(i),forcon(i),fordepth(i)
 992         enddo
 993       endif
 994       write (iout,*)
 995       call hpb_partition
 996       call flush(iout)
 997       return
 998       end