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