source/wham/src-M-SAXS-homology/readrtns.F

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