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