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