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