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