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

   1       subroutine read_control
   2 C
   3 C Read molecular data
   4 C
   5       implicit none
   6       include 'DIMENSIONS'
   7       include 'sizesclu.dat'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.TIME1'
  10       include 'COMMON.SBRIDGE'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.HEADER'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.FREE'
  17       include 'COMMON.INTERACT'
  18       include "COMMON.SPLITELE"
  19       include 'COMMON.SHIELD'
  20       character*320 controlcard,ucase
  21 #ifdef MPL
  22       include 'COMMON.INFO'
  23 #endif
  24       integer i,i1,i2,it1,it2
  25       double precision pi
  26       read (INP,'(a80)') titel
  27       call card_concat(controlcard)
  28
  29       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
  30       call readi(controlcard,'TORMODE',tor_mode,0)
  31       call readi(controlcard,'NRES',nres,0)
  32       call readi(controlcard,'RESCALE',rescale_mode,2)
  33       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
  34       write (iout,*) "DISTCHAINMAX",distchainmax
  35 C Reading the dimensions of box in x,y,z coordinates
  36       call reada(controlcard,'BOXX',boxxsize,100.0d0)
  37       call reada(controlcard,'BOXY',boxysize,100.0d0)
  38       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
  39 c Cutoff range for interactions
  40       call reada(controlcard,"R_CUT",r_cut,15.0d0)
  41       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
  42       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
  43       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
  44       if (lipthick.gt.0.0d0) then
  45        bordliptop=(boxzsize+lipthick)/2.0
  46        bordlipbot=bordliptop-lipthick
  47 C      endif
  48       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
  49      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
  50       buflipbot=bordlipbot+lipbufthick
  51       bufliptop=bordliptop-lipbufthick
  52       if ((lipbufthick*2.0d0).gt.lipthick)
  53      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
  54       endif
  55       write(iout,*) "bordliptop=",bordliptop
  56       write(iout,*) "bordlipbot=",bordlipbot
  57       write(iout,*) "bufliptop=",bufliptop
  58       write(iout,*) "buflipbot=",buflipbot
  59 C Shielding mode
  60       call readi(controlcard,'SHIELD',shield_mode,0)
  61       write (iout,*) "SHIELD MODE",shield_mode
  62       if (shield_mode.gt.0) then
  63       pdbref=(index(controlcard,'PDBREF').gt.0)
  64       if (index(controlcard,"CASC").gt.0) then
  65         iz_sc=1
  66       else if (index(controlcard,"SCONLY").gt.0) then
  67         iz_sc=2
  68       else
  69         iz_sc=0
  70       endif
  71       pi=3.141592d0
  72 C VSolvSphere the volume of solving sphere
  73 C      print *,pi,"pi"
  74 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
  75 C there will be no distinction between proline peptide group and normal peptide
  76 C group in case of shielding parameters
  77       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
  78       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
  79       write (iout,*) VSolvSphere,VSolvSphere_div
  80 C long axis of side chain
  81       do i=1,ntyp
  82       long_r_sidechain(i)=vbldsc0(1,i)
  83       short_r_sidechain(i)=sigma0(i)
  84       enddo
  85       buff_shield=1.0d0
  86       endif
  87       call readi(controlcard,'PDBOUT',outpdb,0)
  88       call readi(controlcard,'MOL2OUT',outmol2,0)
  89       refstr=(index(controlcard,'REFSTR').gt.0)
  90       pdbref=(index(controlcard,'PDBREF').gt.0)
  91       refstr = refstr .or. pdbref
  92       write (iout,*) "REFSTR",refstr," PDBREF",pdbref
  93       iscode=index(controlcard,'ONE_LETTER')
  94       tree=(index(controlcard,'MAKE_TREE').gt.0)
  95       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  96       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  97       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  98      & " CONSTR_DIST",constr_dist
  99       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 100       write (iout,*) "with_theta_constr ",with_theta_constr
 101       call flush(iout)
 102       min_var=(index(controlcard,'MINVAR').gt.0)
 103       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
 104       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
 105       print_fittest=(index(controlcard,'PRINT_FITTEST').gt.0)
 106       call readi(controlcard,'NCUT',ncut,0)
 107       if (ncut.gt.0) then
 108       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
 109       nclust=0
 110       else
 111       call readi(controlcard,'NCLUST',nclust,5)
 112       endif
 113       call readi(controlcard,'SYM',symetr,1)
 114       write (iout,*) 'sym', symetr
 115       call readi(controlcard,'NSTART',nstart,0)
 116       call readi(controlcard,'NEND',nend,0)
 117       call reada(controlcard,'ECUT',ecut,10.0d0)
 118       call reada(controlcard,'PROB',prob_limit,0.99d0)
 119       write (iout,*) "Probability limit",prob_limit
 120       lgrp=(index(controlcard,'LGRP').gt.0)
 121       caonly=(index(controlcard,'CA_ONLY').gt.0)
 122       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
 123       call readi(controlcard,'IOPT',iopt,2)
 124       lefree = index(controlcard,"EFREE").gt.0
 125       call readi(controlcard,'NTEMP',nT,1)
 126       write (iout,*) "nT",nT
 127       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
 128       write (iout,*) "nT",nT
 129       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
 130       do i=1,nT
 131         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
 132       enddo
 133       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
 134       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
 135       lprint_int=index(controlcard,"PRINT_INT") .gt.0
 136       call readi(controlcard,'NSAXS',nsaxs,0)
 137       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 138       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 139       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 140       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 141      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 142       if (min_var) iopt=1
 143       return
 144       end
 145 c--------------------------------------------------------------------------
 146       subroutine molread
 147 C
 148 C Read molecular data.
 149 C
 150       implicit none
 151       include 'DIMENSIONS'
 152       include 'COMMON.IOUNITS'
 153       include 'COMMON.GEO'
 154       include 'COMMON.VAR'
 155       include 'COMMON.INTERACT'
 156       include 'COMMON.LOCAL'
 157       include 'COMMON.NAMES'
 158       include 'COMMON.CHAIN'
 159       include 'COMMON.FFIELD'
 160       include 'COMMON.SBRIDGE'
 161       include 'COMMON.HEADER'
 162       include 'COMMON.CONTROL'
 163       include 'COMMON.CONTACTS'
 164       include 'COMMON.TIME1'
 165       include 'COMMON.TORCNSTR'
 166       include 'COMMON.SHIELD'
 167 #ifdef MPL
 168       include 'COMMON.INFO'
 169 #endif
 170       character*4 sequence(maxres)
 171       character*800 weightcard,controlcard
 172       integer rescode
 173       double precision x(maxvar)
 174       double precision phihel,phibet,sigmahel,sigmabet,sumv,
 175      & secprob(3,maxres)
 176       integer itype_pdb(maxres)
 177       logical seq_comp
 178       integer i,j,kkk,i1,i2,it1,it2,tperm,ii,iperm
 179 C
 180 C Body
 181 C
 182 C Read weights of the subsequent energy terms.
 183       call card_concat(weightcard)
 184       call reada(weightcard,'WSC',wsc,1.0d0)
 185       call reada(weightcard,'WLONG',wsc,wsc)
 186       call reada(weightcard,'WSCP',wscp,1.0d0)
 187       call reada(weightcard,'WELEC',welec,1.0D0)
 188       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 189       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 190       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 191       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 192       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 193       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 194       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 195       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 196       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 197       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 198       call reada(weightcard,'WBOND',wbond,1.0D0)
 199       call reada(weightcard,'WTOR',wtor,1.0D0)
 200       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 201       call reada(weightcard,'WANG',wang,1.0D0)
 202       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 203       call reada(weightcard,'WSAXS',wsaxs,0.0D0)
 204       call reada(weightcard,'SCAL14',scal14,0.4D0)
 205       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 206       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 207       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 208       if (index(weightcard,'SOFT').gt.0) ipot=6
 209       call reada(weightcard,"D0CM",d0cm,3.78d0)
 210       call reada(weightcard,"AKCM",akcm,15.1d0)
 211       call reada(weightcard,"AKTH",akth,11.0d0)
 212       call reada(weightcard,"AKCT",akct,12.0d0)
 213       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 214       call reada(weightcard,"V2SS",v2ss,7.61d0)
 215       call reada(weightcard,"V3SS",v3ss,13.7d0)
 216       call reada(weightcard,"EBR",ebr,-5.50D0)
 217       call reada(weightcard,'WSHIELD',wshield,1.0d0)
 218       write(iout,*) 'WSHIELD',wshield
 219        call reada(weightcard,'WLT',wliptran,0.0D0)
 220       call reada(weightcard,"ATRISS",atriss,0.301D0)
 221       call reada(weightcard,"BTRISS",btriss,0.021D0)
 222       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 223       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 224       write (iout,*) "ATRISS=", atriss
 225       write (iout,*) "BTRISS=", btriss
 226       write (iout,*) "CTRISS=", ctriss
 227       write (iout,*) "DTRISS=", dtriss
 228       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 229       do i=1,maxres
 230         dyn_ss_mask(i)=.false.
 231       enddo
 232       do i=1,maxres-1
 233         do j=i+1,maxres
 234           dyn_ssbond_ij(i,j)=1.0d300
 235         enddo
 236       enddo
 237       call reada(weightcard,"HT",Ht,0.0D0)
 238       if (dyn_ss) then
 239         ss_depth=ebr/wsc-0.25*eps(1,1)
 240         Ht=Ht/wsc-0.25*eps(1,1)
 241         akcm=akcm*wstrain/wsc
 242         akth=akth*wstrain/wsc
 243         akct=akct*wstrain/wsc
 244         v1ss=v1ss*wstrain/wsc
 245         v2ss=v2ss*wstrain/wsc
 246         v3ss=v3ss*wstrain/wsc
 247       else
 248         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 249       endif
 250       write (iout,'(/a)') "Disulfide bridge parameters:"
 251       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 252       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 253       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 254       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 255      & ' v3ss:',v3ss
 256
 257 C 12/1/95 Added weight for the multi-body term WCORR
 258       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 259       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 260       weights(1)=wsc
 261       weights(2)=wscp
 262       weights(3)=welec
 263       weights(4)=wcorr
 264       weights(5)=wcorr5
 265       weights(6)=wcorr6
 266       weights(7)=wel_loc
 267       weights(8)=wturn3
 268       weights(9)=wturn4
 269       weights(10)=wturn6
 270       weights(11)=wang
 271       weights(12)=wscloc
 272       weights(13)=wtor
 273       weights(14)=wtor_d
 274       weights(15)=wstrain
 275       weights(16)=wvdwpp
 276       weights(17)=wbond
 277       weights(18)=scal14
 278       weights(19)=wsccor
 279       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 280      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 281      &  wturn4,wturn6,wsccor
 282    10 format (/'Energy-term weights (unscaled):'//
 283      & 'WSCC=   ',f10.6,' (SC-SC)'/
 284      & 'WSCP=   ',f10.6,' (SC-p)'/
 285      & 'WELEC=  ',f10.6,' (p-p electr)'/
 286      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 287      & 'WBOND=  ',f10.6,' (stretching)'/
 288      & 'WANG=   ',f10.6,' (bending)'/
 289      & 'WSCLOC= ',f10.6,' (SC local)'/
 290      & 'WTOR=   ',f10.6,' (torsional)'/
 291      & 'WTORD=  ',f10.6,' (double torsional)'/
 292      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 293      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 294      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 295      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 296      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 297      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 298      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 299      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 300      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 301
 302       if (wcorr4.gt.0.0d0) then
 303         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 304      &   'between contact pairs of peptide groups'
 305         write (iout,'(2(a,f5.3/))')
 306      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 307      &  'Range of quenching the correlation terms:',2*delt_corr
 308       else if (wcorr.gt.0.0d0) then
 309         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 310      &   'between contact pairs of peptide groups'
 311       endif
 312       write (iout,'(a,f8.3)')
 313      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 314       write (iout,'(a,f8.3)')
 315      &  'General scaling factor of SC-p interactions:',scalscp
 316       r0_corr=cutoff_corr-delt_corr
 317       do i=1,20
 318         aad(i,1)=scalscp*aad(i,1)
 319         aad(i,2)=scalscp*aad(i,2)
 320         bad(i,1)=scalscp*bad(i,1)
 321         bad(i,2)=scalscp*bad(i,2)
 322       enddo
 323
 324       call flush(iout)
 325 c      print *,'indpdb=',indpdb,' pdbref=',pdbref
 326
 327 C Read sequence if not taken from the pdb file.
 328       if (iscode.gt.0) then
 329         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 330       else
 331         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 332       endif
 333 C Convert sequence to numeric code
 334       do i=1,nres
 335         itype(i)=rescode(i,sequence(i),iscode)
 336       enddo
 337 c      print *,nres
 338 c      print '(20i4)',(itype(i),i=1,nres)
 339
 340       do i=1,nres
 341 #ifdef PROCOR
 342         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 343 #else
 344         if (itype(i).eq.ntyp1) then
 345 #endif
 346           itel(i)=0
 347 #ifdef PROCOR
 348         else if (iabs(itype(i+1)).ne.20) then
 349 #else
 350         else if (iabs(itype(i)).ne.20) then
 351 #endif
 352           itel(i)=1
 353         else
 354           itel(i)=2
 355         endif
 356       enddo
 357       write (iout,*) "ITEL"
 358       do i=1,nres-1
 359         write (iout,*) i,itype(i),itel(i)
 360       enddo
 361
 362 c      print *,'Call Read_Bridge.'
 363       call read_bridge
 364 C this fragment reads diheadral constrains
 365       nnt=1
 366       nct=nres
 367 c      print *,'NNT=',NNT,' NCT=',NCT
 368       call seq2chains(nres,itype,nchain,chain_length,chain_border,
 369      &  ireschain)
 370       write(iout,*) "nres",nres," nchain",nchain
 371       do i=1,nchain
 372         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
 373      &    chain_border(2,i)
 374       enddo
 375       call chain_symmetry(nchain,nres,itype,chain_border,
 376      &    chain_length,npermchain,tabpermchain)
 377       do i=1,nres
 378         write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 379      &    ii=1,npermchain)
 380       enddo
 381       write(iout,*) "residue permutations"
 382       do i=1,nres
 383         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
 384       enddo
 385       if (itype(1).eq.ntyp1) nnt=2
 386       if (itype(nres).eq.ntyp1) nct=nct-1
 387       if (nstart.lt.nnt) nstart=nnt
 388       if (nend.gt.nct .or. nend.eq.0) nend=nct
 389       write (iout,*) "nstart",nstart," nend",nend
 390       nres0=nres
 391       if (with_dihed_constr) then
 392
 393       read (inp,*) ndih_constr
 394       if (ndih_constr.gt.0) then
 395         raw_psipred=.false.
 396 C        read (inp,*) ftors
 397 C        write (iout,*) 'FTORS',ftors
 398 C ftors is the force constant for torsional quartic constrains
 399         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 400      &                i=1,ndih_constr)
 401         write (iout,*)
 402      &   'There are',ndih_constr,' constraints on phi angles.'
 403         do i=1,ndih_constr
 404           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 405      &  ftors(i)
 406         enddo
 407         do i=1,ndih_constr
 408           phi0(i)=deg2rad*phi0(i)
 409           drange(i)=deg2rad*drange(i)
 410         enddo
 411       else if (ndih_constr.lt.0) then
 412         raw_psipred=.true.
 413         call card_concat(controlcard)
 414         call reada(controlcard,"PHIHEL",phihel,50.0D0)
 415         call reada(controlcard,"PHIBET",phibet,180.0D0)
 416         call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
 417         call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
 418         call reada(controlcard,"WDIHC",wdihc,0.591d0)
 419         write (iout,*) "Weight of the dihedral restraint term",wdihc
 420         read(inp,'(9x,3f7.3)')
 421      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 422         write (iout,*) "The secprob array"
 423         do i=nnt,nct
 424           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
 425         enddo
 426         ndih_constr=0
 427         do i=nnt+3,nct
 428           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
 429      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
 430             ndih_constr=ndih_constr+1
 431             idih_constr(ndih_constr)=i
 432             sumv=0.0d0
 433             do j=1,3
 434               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
 435               sumv=sumv+vpsipred(j,ndih_constr)
 436             enddo
 437             do j=1,3
 438               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
 439             enddo
 440             phibound(1,ndih_constr)=phihel*deg2rad
 441             phibound(2,ndih_constr)=phibet*deg2rad
 442             sdihed(1,ndih_constr)=sigmahel*deg2rad
 443             sdihed(2,ndih_constr)=sigmabet*deg2rad
 444           endif
 445         enddo
 446         write (iout,*)
 447      &   'There are',ndih_constr,
 448      &   ' bimodal restraints on gamma angles.'
 449         do i=1,ndih_constr
 450           write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
 451      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
 452      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
 453      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
 454      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
 455      &      (vpsipred(j,i),j=1,3)
 456         enddo
 457
 458       endif ! endif ndif_constr.gt.0
 459       endif ! with_dihed_constr
 460       if (with_theta_constr) then
 461 C with_theta_constr is keyword allowing for occurance of theta constrains
 462       read (inp,*) ntheta_constr
 463 C ntheta_constr is the number of theta constrains
 464       if (ntheta_constr.gt.0) then
 465 C        read (inp,*) ftors
 466         read (inp,*) (itheta_constr(i),theta_constr0(i),
 467      &  theta_drange(i),for_thet_constr(i),
 468      &  i=1,ntheta_constr)
 469 C the above code reads from 1 to ntheta_constr
 470 C itheta_constr(i) residue i for which is theta_constr
 471 C theta_constr0 the global minimum value
 472 C theta_drange is range for which there is no energy penalty
 473 C for_thet_constr is the force constant for quartic energy penalty
 474 C E=k*x**4
 475          write (iout,*)
 476      &   'There are',ntheta_constr,' constraints on phi angles.'
 477          do i=1,ntheta_constr
 478           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 479      &    theta_drange(i),
 480      &    for_thet_constr(i)
 481          enddo
 482 C        endif
 483         do i=1,ntheta_constr
 484           theta_constr0(i)=deg2rad*theta_constr0(i)
 485           theta_drange(i)=deg2rad*theta_drange(i)
 486         enddo
 487 C        if(me.eq.king.or..not.out1file)
 488 C     &   write (iout,*) 'FTORS',ftors
 489 C        do i=1,ntheta_constr
 490 C          ii = itheta_constr(i)
 491 C          thetabound(1,ii) = phi0(i)-drange(i)
 492 C          thetabound(2,ii) = phi0(i)+drange(i)
 493 C        enddo
 494       endif ! ntheta_constr.gt.0
 495       endif! with_theta_constr
 496       write (iout,*) "nstart",nstart," nend",nend
 497       write (iout,*) "calling read_saxs_consrtr",nsaxs
 498       if (nsaxs.gt.0) call read_saxs_constr
 499
 500 c      if (pdbref) then
 501 c        read(inp,'(a)') pdbfile
 502 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 503 c        open(ipdbin,file=pdbfile,status='old',err=33)
 504 c        goto 34
 505 c  33    write (iout,'(a)') 'Error opening PDB file.'
 506 c        stop
 507 c  34    continue
 508 c        print *,'Begin reading pdb data'
 509 c        call readpdb
 510 c        print *,'Finished reading pdb data'
 511 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 512 c        do i=1,nres
 513 c          itype_pdb(i)=itype(i)
 514 c        enddo
 515 c        close (ipdbin)
 516 c        write (iout,'(a,i3)') 'nsup=',nsup
 517 c        nstart_seq=nnt
 518 c        if (nsup.le.(nct-nnt+1)) then
 519 c          do i=0,nct-nnt+1-nsup
 520 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 521 c              nstart_seq=nnt+i
 522 c              goto 111
 523 c            endif
 524 c          enddo
 525 c          write (iout,'(a)')
 526 c     &            'Error - sequences to be superposed do not match.'
 527 c          stop
 528 c        else
 529 c          do i=0,nsup-(nct-nnt+1)
 530 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 531 c     &      then
 532 c              nstart_sup=nstart_sup+i
 533 c              nsup=nct-nnt+1
 534 c              goto 111
 535 c            endif
 536 c          enddo
 537 c          write (iout,'(a)')
 538 c     &            'Error - sequences to be superposed do not match.'
 539 c        endif
 540 c  111   continue
 541 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 542 c     &                 ' nstart_seq=',nstart_seq
 543 c      endif
 544       call init_int_table
 545       call setup_var
 546       write (iout,*) "molread: REFSTR",refstr
 547       if (refstr) then
 548         if (.not.pdbref) then
 549           call read_angles(inp,*38)
 550           goto 39
 551    38     write (iout,'(a)') 'Error reading reference structure.'
 552 #ifdef MPL
 553           call mp_stopall(Error_Msg)
 554 #else
 555           stop 'Error reading reference structure'
 556 #endif
 557    39     call chainbuild
 558           nstart_sup=nnt
 559           nstart_seq=nnt
 560           nsup=nct-nnt+1
 561           do i=1,2*nres
 562             do j=1,3
 563               cref(j,i)=c(j,i)
 564             enddo
 565           enddo
 566         endif
 567 c        call contact(.true.,ncont_ref,icont_ref)
 568       endif
 569        if (ns.gt.0) then
 570 C        write (iout,'(/a,i3,a)')
 571 C     &  'The chain contains',ns,' disulfide-bridging cysteines.'
 572         write (iout,'(20i4)') (iss(i),i=1,ns)
 573        if (dyn_ss) then
 574           write(iout,*)"Running with dynamic disulfide-bond formation"
 575        else
 576         write (iout,'(/a/)') 'Pre-formed links are:'
 577         do i=1,nss
 578           i1=ihpb(i)-nres
 579           i2=jhpb(i)-nres
 580           it1=itype(i1)
 581           it2=itype(i2)
 582           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 583      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 584      &    ebr,forcon(i)
 585         enddo
 586         write (iout,'(a)')
 587        endif
 588       endif
 589       if (ns.gt.0.and.dyn_ss) then
 590           do i=nss+1,nhpb
 591             ihpb(i-nss)=ihpb(i)
 592             jhpb(i-nss)=jhpb(i)
 593             forcon(i-nss)=forcon(i)
 594             dhpb(i-nss)=dhpb(i)
 595           enddo
 596           nhpb=nhpb-nss
 597           nss=0
 598           call hpb_partition
 599           do i=1,ns
 600             dyn_ss_mask(iss(i))=.true.
 601           enddo
 602       endif
 603 c Read distance restraints
 604       if (constr_dist.gt.0) then
 605         call read_dist_constr
 606         call hpb_partition
 607       endif
 608       return
 609       end
 610 c-----------------------------------------------------------------------------
 611       logical function seq_comp(itypea,itypeb,length)
 612       implicit none
 613       integer length,itypea(length),itypeb(length)
 614       integer i
 615       do i=1,length
 616         if (itypea(i).ne.itypeb(i)) then
 617           seq_comp=.false.
 618           return
 619         endif
 620       enddo
 621       seq_comp=.true.
 622       return
 623       end
 624 c-----------------------------------------------------------------------------
 625       subroutine read_bridge
 626 C Read information about disulfide bridges.
 627       implicit none
 628       include 'DIMENSIONS'
 629       include 'COMMON.IOUNITS'
 630       include 'COMMON.GEO'
 631       include 'COMMON.VAR'
 632       include 'COMMON.INTERACT'
 633       include 'COMMON.LOCAL'
 634       include 'COMMON.NAMES'
 635       include 'COMMON.CHAIN'
 636       include 'COMMON.FFIELD'
 637       include 'COMMON.SBRIDGE'
 638       include 'COMMON.HEADER'
 639       include 'COMMON.CONTROL'
 640       include 'COMMON.TIME1'
 641 #ifdef MPL
 642       include 'COMMON.INFO'
 643 #endif
 644       integer i,j
 645 C Read bridging residues.
 646       read (inp,*) ns,(iss(i),i=1,ns)
 647 c      print *,'ns=',ns
 648 C Check whether the specified bridging residues are cystines.
 649       do i=1,ns
 650         if (itype(iss(i)).ne.1) then
 651           write (iout,'(2a,i3,a)')
 652      &   'Do you REALLY think that the residue ',
 653      &    restyp(itype(iss(i))),i,
 654      &   ' can form a disulfide bridge?!!!'
 655           write (*,'(2a,i3,a)')
 656      &   'Do you REALLY think that the residue ',
 657      &   restyp(itype(iss(i))),i,
 658      &   ' can form a disulfide bridge?!!!'
 659 #ifdef MPL
 660          call mp_stopall(error_msg)
 661 #else
 662          stop
 663 #endif
 664         endif
 665       enddo
 666 C Read preformed bridges.
 667       if (ns.gt.0) then
 668       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 669       if (nss.gt.0) then
 670         nhpb=nss
 671 C Check if the residues involved in bridges are in the specified list of
 672 C bridging residues.
 673         do i=1,nss
 674           do j=1,i-1
 675             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 676      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 677               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 678      &      ' contains residues present in other pairs.'
 679               write (*,'(a,i3,a)') 'Disulfide pair',i,
 680      &      ' contains residues present in other pairs.'
 681 #ifdef MPL
 682               call mp_stopall(error_msg)
 683 #else
 684               stop
 685 #endif
 686             endif
 687           enddo
 688           do j=1,ns
 689             if (ihpb(i).eq.iss(j)) goto 10
 690           enddo
 691           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 692    10     continue
 693           do j=1,ns
 694             if (jhpb(i).eq.iss(j)) goto 20
 695           enddo
 696           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 697    20     continue
 698 C          dhpb(i)=dbr
 699 C          forcon(i)=fbr
 700         enddo
 701         do i=1,nss
 702           ihpb(i)=ihpb(i)+nres
 703           jhpb(i)=jhpb(i)+nres
 704         enddo
 705       endif
 706       endif
 707       return
 708       end
 709 c----------------------------------------------------------------------------
 710       subroutine read_angles(kanal,*)
 711       implicit none
 712       include 'DIMENSIONS'
 713       include 'COMMON.GEO'
 714       include 'COMMON.VAR'
 715       include 'COMMON.CHAIN'
 716       include 'COMMON.IOUNITS'
 717       integer i,kanal
 718       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 719       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 720       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 721       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 722       do i=1,nres
 723         theta(i)=deg2rad*theta(i)
 724         phi(i)=deg2rad*phi(i)
 725         alph(i)=deg2rad*alph(i)
 726         omeg(i)=deg2rad*omeg(i)
 727       enddo
 728       return
 729    10 return1
 730       end
 731 c----------------------------------------------------------------------------
 732       subroutine reada(rekord,lancuch,wartosc,default)
 733       implicit none
 734       character*(*) rekord,lancuch
 735       double precision wartosc,default
 736       integer ilen,iread
 737       external ilen
 738       iread=index(rekord,lancuch)
 739       if (iread.eq.0) then
 740         wartosc=default
 741         return
 742       endif
 743       iread=iread+ilen(lancuch)+1
 744       read (rekord(iread:),*) wartosc
 745       return
 746       end
 747 c----------------------------------------------------------------------------
 748       subroutine multreada(rekord,lancuch,tablica,dim,default)
 749       implicit none
 750       integer dim,i
 751       double precision tablica(dim),default
 752       character*(*) rekord,lancuch
 753       integer ilen,iread
 754       external ilen
 755       do i=1,dim
 756         tablica(i)=default
 757       enddo
 758       iread=index(rekord,lancuch)
 759       if (iread.eq.0) return
 760       iread=iread+ilen(lancuch)+1
 761       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 762    10 return
 763       end
 764 c----------------------------------------------------------------------------
 765       subroutine readi(rekord,lancuch,wartosc,default)
 766       implicit none
 767       character*(*) rekord,lancuch
 768       integer wartosc,default
 769       integer ilen,iread
 770       external ilen
 771       iread=index(rekord,lancuch)
 772       if (iread.eq.0) then
 773         wartosc=default
 774         return
 775       endif
 776       iread=iread+ilen(lancuch)+1
 777       read (rekord(iread:),*) wartosc
 778       return
 779       end
 780 C----------------------------------------------------------------------
 781       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 782       implicit none
 783       integer dim,i
 784       integer tablica(dim),default
 785       character*(*) rekord,lancuch
 786       character*80 aux
 787       integer ilen,iread
 788       external ilen
 789       do i=1,dim
 790         tablica(i)=default
 791       enddo
 792       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 793       if (iread.eq.0) return
 794       iread=iread+ilen(lancuch)+1
 795       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 796    10 return
 797       end
 798
 799 c----------------------------------------------------------------------------
 800       subroutine card_concat(card)
 801       include 'DIMENSIONS'
 802       include 'COMMON.IOUNITS'
 803       character*(*) card
 804       character*80 karta,ucase
 805       external ilen
 806       read (inp,'(a)') karta
 807       karta=ucase(karta)
 808       card=' '
 809       do while (karta(80:80).eq.'&')
 810         card=card(:ilen(card)+1)//karta(:79)
 811         read (inp,'(a)') karta
 812         karta=ucase(karta)
 813       enddo
 814       card=card(:ilen(card)+1)//karta
 815       return
 816       end
 817 c----------------------------------------------------------------------------
 818       subroutine openunits
 819       implicit none
 820       include 'DIMENSIONS'
 821 #ifdef MPI
 822       include "mpif.h"
 823       character*3 liczba
 824       include "COMMON.MPI"
 825 #endif
 826       include 'COMMON.IOUNITS'
 827       include 'COMMON.CONTROL'
 828       integer lenpre,lenpot,ilen
 829       external ilen
 830       character*16 cformat,cprint
 831       character*16 ucase
 832       integer lenint,lenout
 833       call getenv('INPUT',prefix)
 834       call getenv('OUTPUT',prefout)
 835       call getenv('INTIN',prefintin)
 836       call getenv('COORD',cformat)
 837       call getenv('PRINTCOOR',cprint)
 838       call getenv('SCRATCHDIR',scratchdir)
 839       from_bx=.true.
 840       from_cx=.false.
 841       if (index(ucase(cformat),'CX').gt.0) then
 842         from_cx=.true.
 843         from_bx=.false.
 844       endif
 845       from_cart=.true.
 846       lenpre=ilen(prefix)
 847       lenout=ilen(prefout)
 848       lenint=ilen(prefintin)
 849 C Get the names and open the input files
 850       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 851 #ifdef MPI
 852       write (liczba,'(bz,i3.3)') me
 853       outname=prefout(:lenout)//'_clust.out_'//liczba
 854 #else
 855       outname=prefout(:lenout)//'_clust.out'
 856 #endif
 857       if (from_bx) then
 858         intinname=prefintin(:lenint)//'.bx'
 859       else if (from_cx) then
 860         intinname=prefintin(:lenint)//'.cx'
 861       else
 862         intinname=prefintin(:lenint)//'.int'
 863       endif
 864       rmsname=prefintin(:lenint)//'.rms'
 865       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 866      &       status='unknown')
 867       open (jrms,file=rmsname,status='unknown')
 868       open(iout,file=outname,status='unknown')
 869 C Get parameter filenames and open the parameter files.
 870       call getenv('BONDPAR',bondname)
 871       open (ibond,file=bondname,status='old')
 872       call getenv('THETPAR',thetname)
 873       open (ithep,file=thetname,status='old')
 874       call getenv('ROTPAR',rotname)
 875       open (irotam,file=rotname,status='old')
 876       call getenv('TORPAR',torname)
 877       open (itorp,file=torname,status='old')
 878 #ifndef NEWCORR
 879       call getenv('TORDPAR',tordname)
 880       open (itordp,file=tordname,status='old')
 881 #endif
 882       call getenv('FOURIER',fouriername)
 883       open (ifourier,file=fouriername,status='old')
 884       call getenv('ELEPAR',elename)
 885       open (ielep,file=elename,status='old')
 886       call getenv('SIDEPAR',sidename)
 887       open (isidep,file=sidename,status='old')
 888       call getenv('SIDEP',sidepname)
 889       open (isidep1,file=sidepname,status="old")
 890       call getenv('SCCORPAR',sccorname)
 891       open (isccor,file=sccorname,status="old")
 892       call getenv('LIPTRANPAR',liptranname)
 893       open (iliptranpar,file=liptranname,status='old')
 894 #ifndef OLDSCP
 895 C
 896 C 8/9/01 In the newest version SCp interaction constants are read from a file
 897 C Use -DOLDSCP to use hard-coded constants instead.
 898 C
 899       call getenv('SCPPAR',scpname)
 900       open (iscpp,file=scpname,status='old')
 901 #endif
 902       return
 903       end
 904 c--------------------------------------------------------------------------
 905       subroutine read_dist_constr
 906       implicit real*8 (a-h,o-z)
 907       include 'DIMENSIONS'
 908 #ifdef MPI
 909       include 'mpif.h'
 910 #endif
 911       include 'COMMON.CONTROL'
 912       include 'COMMON.CHAIN'
 913       include 'COMMON.IOUNITS'
 914       include 'COMMON.SBRIDGE'
 915       include 'COMMON.INTERACT'
 916       integer ifrag_(2,100),ipair_(2,100)
 917       double precision wfrag_(100),wpair_(100)
 918       character*500 controlcard
 919       logical lprn /.true./
 920       logical normalize,next
 921       integer restr_type
 922       double precision scal_bfac
 923       double precision xlink(4,0:4) /
 924 c           a          b       c     sigma
 925      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
 926      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
 927      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
 928      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
 929      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
 930       write (iout,*) "Calling read_dist_constr"
 931 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 932 c      call flush(iout)
 933       next=.true.
 934
 935       DO WHILE (next)
 936
 937       call card_concat(controlcard)
 938       next = index(controlcard,"NEXT").gt.0
 939       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
 940       write (iout,*) "restr_type",restr_type
 941       call readi(controlcard,"NFRAG",nfrag_,0)
 942       call readi(controlcard,"NFRAG",nfrag_,0)
 943       call readi(controlcard,"NPAIR",npair_,0)
 944       call readi(controlcard,"NDIST",ndist_,0)
 945       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 946       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
 947       if (restr_type.eq.10)
 948      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
 949       if (restr_type.eq.12)
 950      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
 951       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 952       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 953       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 954       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 955       normalize = index(controlcard,"NORMALIZE").gt.0
 956       write (iout,*) "WBOLTZD",wboltzd
 957       write (iout,*) "SCAL_PEAK",scal_peak
 958       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 959       write (iout,*) "IFRAG"
 960       do i=1,nfrag_
 961         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 962       enddo
 963       write (iout,*) "IPAIR"
 964       do i=1,npair_
 965         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 966       enddo
 967       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
 968         nres0=nres
 969         read(inp,'(a)') pdbfile
 970         write (iout,*)
 971      & "Distance restraints will be constructed from structure ",pdbfile
 972         open(ipdbin,file=pdbfile,status='old',err=11)
 973         call readpdb(.true.)
 974         nres=nres0
 975         close(ipdbin)
 976       endif
 977       do i=1,nfrag_
 978         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 979         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 980      &    ifrag_(2,i)=nstart_sup+nsup-1
 981 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 982 c        call flush(iout)
 983         if (wfrag_(i).eq.0.0d0) cycle
 984         do j=ifrag_(1,i),ifrag_(2,i)-1
 985           do k=j+1,ifrag_(2,i)
 986 c            write (iout,*) "j",j," k",k
 987             ddjk=dist(j,k)
 988             if (restr_type.eq.1) then
 989               nhpb=nhpb+1
 990               irestr_type(nhpb)=1
 991               ihpb(nhpb)=j
 992               jhpb(nhpb)=k
 993               dhpb(nhpb)=ddjk
 994               forcon(nhpb)=wfrag_(i)
 995             else if (constr_dist.eq.2) then
 996               if (ddjk.le.dist_cut) then
 997                 nhpb=nhpb+1
 998                 irestr_type(nhpb)=1
 999                 ihpb(nhpb)=j
1000                 jhpb(nhpb)=k
1001                 dhpb(nhpb)=ddjk
1002                 forcon(nhpb)=wfrag_(i)
1003               endif
1004             else if (restr_type.eq.3) then
1005               nhpb=nhpb+1
1006               irestr_type(nhpb)=1
1007               ihpb(nhpb)=j
1008               jhpb(nhpb)=k
1009               dhpb(nhpb)=ddjk
1010               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1011             endif
1012             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
1013      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1014           enddo
1015         enddo
1016       enddo
1017       do i=1,npair_
1018         if (wpair_(i).eq.0.0d0) cycle
1019         ii = ipair_(1,i)
1020         jj = ipair_(2,i)
1021         if (ii.gt.jj) then
1022           itemp=ii
1023           ii=jj
1024           jj=itemp
1025         endif
1026         do j=ifrag_(1,ii),ifrag_(2,ii)
1027           do k=ifrag_(1,jj),ifrag_(2,jj)
1028             ddjk=dist(j,k)
1029             if (restr_type.eq.1) then
1030               nhpb=nhpb+1
1031               irestr_type(nhpb)=1
1032               ihpb(nhpb)=j
1033               jhpb(nhpb)=k
1034               dhpb(nhpb)=ddjk
1035               forcon(nhpb)=wpair_(i)
1036             else if (constr_dist.eq.2) then
1037               if (ddjk.le.dist_cut) then
1038                 nhpb=nhpb+1
1039                 irestr_type(nhpb)=1
1040                 ihpb(nhpb)=j
1041                 jhpb(nhpb)=k
1042                 dhpb(nhpb)=ddjk
1043                 forcon(nhpb)=wpair_(i)
1044               endif
1045             else if (restr_type.eq.3) then
1046               nhpb=nhpb+1
1047               irestr_type(nhpb)=1
1048               ihpb(nhpb)=j
1049               jhpb(nhpb)=k
1050               dhpb(nhpb)=ddjk
1051               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1052             endif
1053             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1054      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1055           enddo
1056         enddo
1057       enddo
1058
1059 c      print *,ndist_
1060       write (iout,*) "Distance restraints as read from input"
1061       do i=1,ndist_
1062         if (restr_type.eq.12) then
1063           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1064      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1065      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1066      &    fordepth_peak(nhpb_peak+1),npeak
1067 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1068 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1069 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1070 c     &    fordepth_peak(nhpb_peak+1),npeak
1071           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1072      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1073           nhpb_peak=nhpb_peak+1
1074           irestr_type_peak(nhpb_peak)=12
1075           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1076           ipeak(2,npeak)=i
1077           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1078      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1079      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1080      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1081      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1082           if (ibecarb_peak(nhpb_peak).eq.3) then
1083             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1084           else if (ibecarb_peak(nhpb_peak).eq.2) then
1085             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1086           else if (ibecarb_peak(nhpb_peak).eq.1) then
1087             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1088             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1089           endif
1090         else if (restr_type.eq.11) then
1091           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1092      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1093 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1094           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1095           nhpb=nhpb+1
1096           irestr_type(nhpb)=11
1097           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1098      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1099      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1100 c          if (ibecarb(nhpb).gt.0) then
1101 c            ihpb(nhpb)=ihpb(nhpb)+nres
1102 c            jhpb(nhpb)=jhpb(nhpb)+nres
1103 c          endif
1104          if (ibecarb(nhpb).eq.3) then
1105             ihpb(nhpb)=ihpb(nhpb)+nres
1106           else if (ibecarb(nhpb).eq.2) then
1107             ihpb(nhpb)=ihpb(nhpb)+nres
1108           else if (ibecarb(nhpb).eq.1) then
1109             ihpb(nhpb)=ihpb(nhpb)+nres
1110             jhpb(nhpb)=jhpb(nhpb)+nres
1111           endif
1112         else if (restr_type.eq.10) then
1113 c Cross-lonk Markov-like potential
1114           call card_concat(controlcard)
1115           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1116           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1117           ibecarb(nhpb+1)=0
1118           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1119           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1120           if (index(controlcard,"ZL").gt.0) then
1121             link_type=1
1122           else if (index(controlcard,"ADH").gt.0) then
1123             link_type=2
1124           else if (index(controlcard,"PDH").gt.0) then
1125             link_type=3
1126           else if (index(controlcard,"DSS").gt.0) then
1127             link_type=4
1128           else
1129             link_type=0
1130           endif
1131           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1132      &       xlink(1,link_type))
1133           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1134      &       xlink(2,link_type))
1135           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1136      &       xlink(3,link_type))
1137           call reada(controlcard,"SIGMA",forcon(nhpb+1),
1138      &       xlink(4,link_type))
1139           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1140 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1141 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1142           if (forcon(nhpb+1).le.0.0d0 .or.
1143      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1144           nhpb=nhpb+1
1145           irestr_type(nhpb)=10
1146 c          if (ibecarb(nhpb).gt.0) then
1147 c            ihpb(nhpb)=ihpb(nhpb)+nres
1148 c            jhpb(nhpb)=jhpb(nhpb)+nres
1149 c          endif
1150           if (ibecarb(nhpb).eq.3) then
1151             jhpb(nhpb)=jhpb(nhpb)+nres
1152           else if (ibecarb(nhpb).eq.2) then
1153             ihpb(nhpb)=ihpb(nhpb)+nres
1154           else if (ibecarb(nhpb).eq.1) then
1155             ihpb(nhpb)=ihpb(nhpb)+nres
1156             jhpb(nhpb)=jhpb(nhpb)+nres
1157           endif
1158           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1159      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1160      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1161      &     irestr_type(nhpb)
1162         else
1163 C        print *,"in else"
1164           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1165      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1166           if (forcon(nhpb+1).gt.0.0d0) then
1167           nhpb=nhpb+1
1168           if (dhpb1(nhpb).eq.0.0d0) then
1169             irestr_type(nhpb)=1
1170           else
1171             irestr_type(nhpb)=2
1172           endif
1173 c          if (ibecarb(nhpb).gt.0) then
1174 c            ihpb(nhpb)=ihpb(nhpb)+nres
1175 c            jhpb(nhpb)=jhpb(nhpb)+nres
1176 c          endif
1177           if (ibecarb(nhpb).eq.3) then
1178             jhpb(nhpb)=jhpb(nhpb)+nres
1179           else if (ibecarb(nhpb).eq.2) then
1180             ihpb(nhpb)=ihpb(nhpb)+nres
1181           else if (ibecarb(nhpb).eq.1) then
1182             ihpb(nhpb)=ihpb(nhpb)+nres
1183             jhpb(nhpb)=jhpb(nhpb)+nres
1184           endif
1185           if (dhpb(nhpb).eq.0.0d0)
1186      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1187         endif
1188         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1189      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1190         endif
1191 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1192 C        if (forcon(nhpb+1).gt.0.0d0) then
1193 C          nhpb=nhpb+1
1194 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1195       enddo
1196
1197       if (restr_type.eq.4) then
1198         write (iout,*) "The BFAC array"
1199         do i=nnt,nct
1200           write (iout,'(i5,f10.5)') i,bfac(i)
1201         enddo
1202         do i=nnt,nct
1203           if (itype(i).eq.ntyp1) cycle
1204           do j=nnt,i-1
1205             if (itype(j).eq.ntyp1) cycle
1206             if (itype(i).eq.10) then
1207               iiend=0
1208             else
1209               iiend=1
1210             endif
1211             if (itype(j).eq.10) then
1212               jjend=0
1213             else
1214               jjend=1
1215             endif
1216             kk=0
1217             do ii=0,iiend
1218             do jj=0,jjend
1219             nhpb=nhpb+1
1220             irestr_type(nhpb)=1
1221             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1222             irestr_type(nhpb)=1
1223             ibecarb(nhpb)=kk
1224             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1225             ihpb(nhpb)=i+nres*ii
1226             jhpb(nhpb)=j+nres*jj
1227             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1228             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1229      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1230      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1231      &       irestr_type(nhpb)
1232             kk=kk+1
1233           enddo
1234           enddo
1235           enddo
1236         enddo
1237       endif
1238
1239       if (restr_type.eq.5) then
1240         restr_on_coord=.true.
1241         do i=nnt,nct
1242           if (itype(i).eq.ntyp1) cycle
1243           bfac(i)=(scal_bfac/bfac(i))**2
1244         enddo
1245       endif
1246
1247       ENDDO ! next
1248
1249       fordepthmax=0.0d0
1250       if (normalize) then
1251         do i=nss+1,nhpb
1252           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1253      &      fordepthmax=fordepth(i)
1254         enddo
1255         do i=nss+1,nhpb
1256           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1257         enddo
1258       endif
1259       if (nhpb.gt.nss)  then
1260         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1261      &  "The following",nhpb-nss,
1262      &  " distance restraints have been imposed:",
1263      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
1264      &  "  score"," type"
1265         do i=nss+1,nhpb
1266           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1267      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1268      &  irestr_type(i)
1269         enddo
1270       endif
1271       call hpb_partition
1272       call flush(iout)
1273       return
1274    11 write (iout,*)"read_dist_restr: error reading reference structure"
1275       stop
1276       end
1277 c-------------------------------------------------------------------------------
1278       subroutine read_saxs_constr
1279       implicit real*8 (a-h,o-z)
1280       include 'DIMENSIONS'
1281 #ifdef MPI
1282       include 'mpif.h'
1283 #endif
1284       include 'COMMON.CONTROL'
1285       include 'COMMON.CHAIN'
1286       include 'COMMON.IOUNITS'
1287       include 'COMMON.SBRIDGE'
1288       double precision cm(3)
1289 c      read(inp,*) nsaxs
1290       write (iout,*) "Calling read_saxs nsaxs",nsaxs
1291       call flush(iout)
1292       if (saxs_mode.eq.0) then
1293 c SAXS distance distribution
1294       do i=1,nsaxs
1295         read(inp,*) distsaxs(i),Psaxs(i)
1296       enddo
1297       Cnorm = 0.0d0
1298       do i=1,nsaxs
1299         Cnorm = Cnorm + Psaxs(i)
1300       enddo
1301       write (iout,*) "Cnorm",Cnorm
1302       do i=1,nsaxs
1303         Psaxs(i)=Psaxs(i)/Cnorm
1304       enddo
1305       write (iout,*) "Normalized distance distribution from SAXS"
1306       do i=1,nsaxs
1307         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1308       enddo
1309       Wsaxs0=0.0d0
1310       do i=1,nsaxs
1311         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1312       enddo
1313       write (iout,*) "Wsaxs0",Wsaxs0
1314       else
1315 c SAXS "spheres".
1316       do i=1,nsaxs
1317         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1318       enddo
1319       do j=1,3
1320         cm(j)=0.0d0
1321       enddo
1322       do i=1,nsaxs
1323         do j=1,3
1324           cm(j)=cm(j)+Csaxs(j,i)
1325         enddo
1326       enddo
1327       do j=1,3
1328         cm(j)=cm(j)/nsaxs
1329       enddo
1330       do i=1,nsaxs
1331         do j=1,3
1332           Csaxs(j,i)=Csaxs(j,i)-cm(j)
1333         enddo
1334       enddo
1335       write (iout,*) "SAXS sphere coordinates"
1336       do i=1,nsaxs
1337         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
1338       enddo
1339       endif
1340       return
1341       end