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