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

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