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