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