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,"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,max_cyst-1
 246         do j=i+1,max_cyst
 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 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
 720 c     numbers
 721       icys=0
 722       do i=1,ns
 723         icys(iss(i))=i
 724       enddo
 725 C Check whether the specified bridging residues are cystines.
 726       do i=1,ns
 727         if (itype(iss(i)).ne.1) then
 728           write (iout,'(2a,i3,a)')
 729      &   'Do you REALLY think that the residue ',
 730      &    restyp(itype(iss(i))),i,
 731      &   ' can form a disulfide bridge?!!!'
 732           write (*,'(2a,i3,a)')
 733      &   'Do you REALLY think that the residue ',
 734      &   restyp(itype(iss(i))),i,
 735      &   ' can form a disulfide bridge?!!!'
 736 #ifdef MPL
 737          call mp_stopall(error_msg)
 738 #else
 739          stop
 740 #endif
 741         endif
 742       enddo
 743 C Read preformed bridges.
 744       if (ns.gt.0) then
 745       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 746       if (nss.gt.0) then
 747         nhpb=nss
 748 C Check if the residues involved in bridges are in the specified list of
 749 C bridging residues.
 750         do i=1,nss
 751           do j=1,i-1
 752             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 753      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 754               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 755      &      ' contains residues present in other pairs.'
 756               write (*,'(a,i3,a)') 'Disulfide pair',i,
 757      &      ' contains residues present in other pairs.'
 758 #ifdef MPL
 759               call mp_stopall(error_msg)
 760 #else
 761               stop
 762 #endif
 763             endif
 764           enddo
 765           do j=1,ns
 766             if (ihpb(i).eq.iss(j)) goto 10
 767           enddo
 768           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 769    10     continue
 770           do j=1,ns
 771             if (jhpb(i).eq.iss(j)) goto 20
 772           enddo
 773           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 774    20     continue
 775 C          dhpb(i)=dbr
 776 C          forcon(i)=fbr
 777         enddo
 778         do i=1,nss
 779           ihpb(i)=ihpb(i)+nres
 780           jhpb(i)=jhpb(i)+nres
 781         enddo
 782       endif
 783       endif
 784       return
 785       end
 786 c----------------------------------------------------------------------------
 787       subroutine read_angles(kanal,*)
 788       implicit none
 789       include 'DIMENSIONS'
 790       include 'COMMON.GEO'
 791       include 'COMMON.VAR'
 792       include 'COMMON.CHAIN'
 793       include 'COMMON.IOUNITS'
 794       integer i,kanal
 795       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 796       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 797       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 798       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 799       do i=1,nres
 800         theta(i)=deg2rad*theta(i)
 801         phi(i)=deg2rad*phi(i)
 802         alph(i)=deg2rad*alph(i)
 803         omeg(i)=deg2rad*omeg(i)
 804       enddo
 805       return
 806    10 return1
 807       end
 808 c----------------------------------------------------------------------------
 809       subroutine reada(rekord,lancuch,wartosc,default)
 810       implicit none
 811       character*(*) rekord,lancuch
 812       double precision wartosc,default
 813       integer ilen,iread
 814       external ilen
 815       iread=index(rekord,lancuch)
 816       if (iread.eq.0) then
 817         wartosc=default
 818         return
 819       endif
 820       iread=iread+ilen(lancuch)+1
 821       read (rekord(iread:),*) wartosc
 822       return
 823       end
 824 c----------------------------------------------------------------------------
 825       subroutine multreada(rekord,lancuch,tablica,dim,default)
 826       implicit none
 827       integer dim,i
 828       double precision tablica(dim),default
 829       character*(*) rekord,lancuch
 830       integer ilen,iread
 831       external ilen
 832       do i=1,dim
 833         tablica(i)=default
 834       enddo
 835       iread=index(rekord,lancuch)
 836       if (iread.eq.0) return
 837       iread=iread+ilen(lancuch)+1
 838       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 839    10 return
 840       end
 841 c----------------------------------------------------------------------------
 842       subroutine readi(rekord,lancuch,wartosc,default)
 843       implicit none
 844       character*(*) rekord,lancuch
 845       integer wartosc,default
 846       integer ilen,iread
 847       external ilen
 848       iread=index(rekord,lancuch)
 849       if (iread.eq.0) then
 850         wartosc=default
 851         return
 852       endif
 853       iread=iread+ilen(lancuch)+1
 854       read (rekord(iread:),*) wartosc
 855       return
 856       end
 857 C----------------------------------------------------------------------
 858       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 859       implicit none
 860       integer dim,i
 861       integer tablica(dim),default
 862       character*(*) rekord,lancuch
 863       character*80 aux
 864       integer ilen,iread
 865       external ilen
 866       do i=1,dim
 867         tablica(i)=default
 868       enddo
 869       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 870       if (iread.eq.0) return
 871       iread=iread+ilen(lancuch)+1
 872       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 873    10 return
 874       end
 875
 876 c----------------------------------------------------------------------------
 877       subroutine card_concat(card)
 878       include 'DIMENSIONS'
 879       include 'COMMON.IOUNITS'
 880       character*(*) card
 881       character*80 karta,ucase
 882       external ilen
 883       read (inp,'(a)') karta
 884       karta=ucase(karta)
 885       card=' '
 886       do while (karta(80:80).eq.'&')
 887         card=card(:ilen(card)+1)//karta(:79)
 888         read (inp,'(a)') karta
 889         karta=ucase(karta)
 890       enddo
 891       card=card(:ilen(card)+1)//karta
 892       return
 893       end
 894 c----------------------------------------------------------------------------
 895       subroutine openunits
 896       implicit none
 897       include 'DIMENSIONS'
 898 #ifdef MPI
 899       include "mpif.h"
 900       character*3 liczba
 901       include "COMMON.MPI"
 902 #endif
 903       include 'COMMON.IOUNITS'
 904       include 'COMMON.CONTROL'
 905       integer lenpre,lenpot,ilen
 906       external ilen
 907       character*16 cformat,cprint
 908       character*16 ucase
 909       integer lenint,lenout
 910       call getenv('INPUT',prefix)
 911       call getenv('OUTPUT',prefout)
 912       call getenv('INTIN',prefintin)
 913       call getenv('COORD',cformat)
 914       call getenv('PRINTCOOR',cprint)
 915       call getenv('SCRATCHDIR',scratchdir)
 916       from_bx=.true.
 917       from_cx=.false.
 918       if (index(ucase(cformat),'CX').gt.0) then
 919         from_cx=.true.
 920         from_bx=.false.
 921       endif
 922       from_cart=.true.
 923       lenpre=ilen(prefix)
 924       lenout=ilen(prefout)
 925       lenint=ilen(prefintin)
 926 C Get the names and open the input files
 927       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 928 #ifdef MPI
 929       write (liczba,'(bz,i3.3)') me
 930       outname=prefout(:lenout)//'_clust.out_'//liczba
 931 #else
 932       outname=prefout(:lenout)//'_clust.out'
 933 #endif
 934       if (from_bx) then
 935         intinname=prefintin(:lenint)//'.bx'
 936       else if (from_cx) then
 937         intinname=prefintin(:lenint)//'.cx'
 938       else
 939         intinname=prefintin(:lenint)//'.int'
 940       endif
 941       rmsname=prefintin(:lenint)//'.rms'
 942       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 943      &       status='unknown')
 944       open (jrms,file=rmsname,status='unknown')
 945       open(iout,file=outname,status='unknown')
 946 C Get parameter filenames and open the parameter files.
 947       call getenv('BONDPAR',bondname)
 948       open (ibond,file=bondname,status='old')
 949       call getenv('THETPAR',thetname)
 950       open (ithep,file=thetname,status='old')
 951       call getenv('ROTPAR',rotname)
 952       open (irotam,file=rotname,status='old')
 953       call getenv('TORPAR',torname)
 954       open (itorp,file=torname,status='old')
 955 #ifndef NEWCORR
 956       call getenv('TORDPAR',tordname)
 957       open (itordp,file=tordname,status='old')
 958 #endif
 959       call getenv('FOURIER',fouriername)
 960       open (ifourier,file=fouriername,status='old')
 961       call getenv('ELEPAR',elename)
 962       open (ielep,file=elename,status='old')
 963       call getenv('SIDEPAR',sidename)
 964       open (isidep,file=sidename,status='old')
 965       call getenv('SIDEP',sidepname)
 966       open (isidep1,file=sidepname,status="old")
 967       call getenv('SCCORPAR',sccorname)
 968       open (isccor,file=sccorname,status="old")
 969       call getenv('LIPTRANPAR',liptranname)
 970       open (iliptranpar,file=liptranname,status='old')
 971 #ifndef OLDSCP
 972 C
 973 C 8/9/01 In the newest version SCp interaction constants are read from a file
 974 C Use -DOLDSCP to use hard-coded constants instead.
 975 C
 976       call getenv('SCPPAR',scpname)
 977       open (iscpp,file=scpname,status='old')
 978 #endif
 979       return
 980       end
 981 c--------------------------------------------------------------------------
 982       subroutine read_dist_constr
 983       implicit real*8 (a-h,o-z)
 984       include 'DIMENSIONS'
 985 #ifdef MPI
 986       include 'mpif.h'
 987 #endif
 988       include 'COMMON.CONTROL'
 989       include 'COMMON.CHAIN'
 990       include 'COMMON.IOUNITS'
 991       include 'COMMON.SBRIDGE'
 992       include 'COMMON.INTERACT'
 993       integer ifrag_(2,100),ipair_(2,100)
 994       double precision wfrag_(100),wpair_(100)
 995       character*500 controlcard
 996       logical lprn /.true./
 997       logical normalize,next
 998       integer restr_type
 999       double precision scal_bfac
1000       double precision xlink(4,0:4) /
1001 c           a          b       c     sigma
1002      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
1003      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
1004      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
1005      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
1006      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
1007       write (iout,*) "Calling read_dist_constr"
1008 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1009 c      call flush(iout)
1010       next=.true.
1011
1012       DO WHILE (next)
1013
1014       call card_concat(controlcard)
1015       next = index(controlcard,"NEXT").gt.0
1016       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
1017       write (iout,*) "restr_type",restr_type
1018       call readi(controlcard,"NFRAG",nfrag_,0)
1019       call readi(controlcard,"NFRAG",nfrag_,0)
1020       call readi(controlcard,"NPAIR",npair_,0)
1021       call readi(controlcard,"NDIST",ndist_,0)
1022       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1023       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
1024       if (restr_type.eq.10)
1025      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
1026       if (restr_type.eq.12)
1027      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
1028       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1029       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1030       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1031       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1032       normalize = index(controlcard,"NORMALIZE").gt.0
1033       write (iout,*) "WBOLTZD",wboltzd
1034       write (iout,*) "SCAL_PEAK",scal_peak
1035       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1036       write (iout,*) "IFRAG"
1037       do i=1,nfrag_
1038         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1039       enddo
1040       write (iout,*) "IPAIR"
1041       do i=1,npair_
1042         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1043       enddo
1044       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
1045         nres0=nres
1046         read(inp,'(a)') pdbfile
1047         write (iout,*)
1048      & "Distance restraints will be constructed from structure ",pdbfile
1049         open(ipdbin,file=pdbfile,status='old',err=11)
1050         call readpdb(.true.)
1051         nres=nres0
1052         close(ipdbin)
1053       endif
1054       do i=1,nfrag_
1055         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1056         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1057      &    ifrag_(2,i)=nstart_sup+nsup-1
1058 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1059 c        call flush(iout)
1060         if (wfrag_(i).eq.0.0d0) cycle
1061         do j=ifrag_(1,i),ifrag_(2,i)-1
1062           do k=j+1,ifrag_(2,i)
1063 c            write (iout,*) "j",j," k",k
1064             ddjk=dist(j,k)
1065             if (restr_type.eq.1) then
1066               nhpb=nhpb+1
1067               irestr_type(nhpb)=1
1068               ihpb(nhpb)=j
1069               jhpb(nhpb)=k
1070               dhpb(nhpb)=ddjk
1071               forcon(nhpb)=wfrag_(i)
1072             else if (constr_dist.eq.2) then
1073               if (ddjk.le.dist_cut) then
1074                 nhpb=nhpb+1
1075                 irestr_type(nhpb)=1
1076                 ihpb(nhpb)=j
1077                 jhpb(nhpb)=k
1078                 dhpb(nhpb)=ddjk
1079                 forcon(nhpb)=wfrag_(i)
1080               endif
1081             else if (restr_type.eq.3) then
1082               nhpb=nhpb+1
1083               irestr_type(nhpb)=1
1084               ihpb(nhpb)=j
1085               jhpb(nhpb)=k
1086               dhpb(nhpb)=ddjk
1087               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1088             endif
1089             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
1090      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1091           enddo
1092         enddo
1093       enddo
1094       do i=1,npair_
1095         if (wpair_(i).eq.0.0d0) cycle
1096         ii = ipair_(1,i)
1097         jj = ipair_(2,i)
1098         if (ii.gt.jj) then
1099           itemp=ii
1100           ii=jj
1101           jj=itemp
1102         endif
1103         do j=ifrag_(1,ii),ifrag_(2,ii)
1104           do k=ifrag_(1,jj),ifrag_(2,jj)
1105             ddjk=dist(j,k)
1106             if (restr_type.eq.1) then
1107               nhpb=nhpb+1
1108               irestr_type(nhpb)=1
1109               ihpb(nhpb)=j
1110               jhpb(nhpb)=k
1111               dhpb(nhpb)=ddjk
1112               forcon(nhpb)=wpair_(i)
1113             else if (constr_dist.eq.2) then
1114               if (ddjk.le.dist_cut) then
1115                 nhpb=nhpb+1
1116                 irestr_type(nhpb)=1
1117                 ihpb(nhpb)=j
1118                 jhpb(nhpb)=k
1119                 dhpb(nhpb)=ddjk
1120                 forcon(nhpb)=wpair_(i)
1121               endif
1122             else if (restr_type.eq.3) then
1123               nhpb=nhpb+1
1124               irestr_type(nhpb)=1
1125               ihpb(nhpb)=j
1126               jhpb(nhpb)=k
1127               dhpb(nhpb)=ddjk
1128               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1129             endif
1130             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1131      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1132           enddo
1133         enddo
1134       enddo
1135
1136 c      print *,ndist_
1137       write (iout,*) "Distance restraints as read from input"
1138       do i=1,ndist_
1139         if (restr_type.eq.12) then
1140           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1141      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1142      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1143      &    fordepth_peak(nhpb_peak+1),npeak
1144 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1145 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1146 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1147 c     &    fordepth_peak(nhpb_peak+1),npeak
1148           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1149      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1150           nhpb_peak=nhpb_peak+1
1151           irestr_type_peak(nhpb_peak)=12
1152           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1153           ipeak(2,npeak)=i
1154           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1155      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1156      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1157      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1158      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1159           if (ibecarb_peak(nhpb_peak).eq.3) then
1160             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1161           else if (ibecarb_peak(nhpb_peak).eq.2) then
1162             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1163           else if (ibecarb_peak(nhpb_peak).eq.1) then
1164             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1165             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1166           endif
1167         else if (restr_type.eq.11) then
1168           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1169      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1170 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1171           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1172           nhpb=nhpb+1
1173           irestr_type(nhpb)=11
1174           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1175      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1176      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1177 c          if (ibecarb(nhpb).gt.0) then
1178 c            ihpb(nhpb)=ihpb(nhpb)+nres
1179 c            jhpb(nhpb)=jhpb(nhpb)+nres
1180 c          endif
1181          if (ibecarb(nhpb).eq.3) then
1182             ihpb(nhpb)=ihpb(nhpb)+nres
1183           else if (ibecarb(nhpb).eq.2) then
1184             ihpb(nhpb)=ihpb(nhpb)+nres
1185           else if (ibecarb(nhpb).eq.1) then
1186             ihpb(nhpb)=ihpb(nhpb)+nres
1187             jhpb(nhpb)=jhpb(nhpb)+nres
1188           endif
1189         else if (restr_type.eq.10) then
1190 c Cross-lonk Markov-like potential
1191           call card_concat(controlcard)
1192           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1193           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1194           ibecarb(nhpb+1)=0
1195           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1196           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1197           if (index(controlcard,"ZL").gt.0) then
1198             link_type=1
1199           else if (index(controlcard,"ADH").gt.0) then
1200             link_type=2
1201           else if (index(controlcard,"PDH").gt.0) then
1202             link_type=3
1203           else if (index(controlcard,"DSS").gt.0) then
1204             link_type=4
1205           else
1206             link_type=0
1207           endif
1208           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1209      &       xlink(1,link_type))
1210           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1211      &       xlink(2,link_type))
1212           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1213      &       xlink(3,link_type))
1214           call reada(controlcard,"SIGMA",forcon(nhpb+1),
1215      &       xlink(4,link_type))
1216           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1217 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1218 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1219           if (forcon(nhpb+1).le.0.0d0 .or.
1220      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1221           nhpb=nhpb+1
1222           irestr_type(nhpb)=10
1223 c          if (ibecarb(nhpb).gt.0) then
1224 c            ihpb(nhpb)=ihpb(nhpb)+nres
1225 c            jhpb(nhpb)=jhpb(nhpb)+nres
1226 c          endif
1227           if (ibecarb(nhpb).eq.3) then
1228             jhpb(nhpb)=jhpb(nhpb)+nres
1229           else if (ibecarb(nhpb).eq.2) then
1230             ihpb(nhpb)=ihpb(nhpb)+nres
1231           else if (ibecarb(nhpb).eq.1) then
1232             ihpb(nhpb)=ihpb(nhpb)+nres
1233             jhpb(nhpb)=jhpb(nhpb)+nres
1234           endif
1235           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1236      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1237      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1238      &     irestr_type(nhpb)
1239         else
1240 C        print *,"in else"
1241           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1242      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1243           if (forcon(nhpb+1).gt.0.0d0) then
1244           nhpb=nhpb+1
1245           if (dhpb1(nhpb).eq.0.0d0) then
1246             irestr_type(nhpb)=1
1247           else
1248             irestr_type(nhpb)=2
1249           endif
1250 c          if (ibecarb(nhpb).gt.0) then
1251 c            ihpb(nhpb)=ihpb(nhpb)+nres
1252 c            jhpb(nhpb)=jhpb(nhpb)+nres
1253 c          endif
1254           if (ibecarb(nhpb).eq.3) then
1255             jhpb(nhpb)=jhpb(nhpb)+nres
1256           else if (ibecarb(nhpb).eq.2) then
1257             ihpb(nhpb)=ihpb(nhpb)+nres
1258           else if (ibecarb(nhpb).eq.1) then
1259             ihpb(nhpb)=ihpb(nhpb)+nres
1260             jhpb(nhpb)=jhpb(nhpb)+nres
1261           endif
1262           if (dhpb(nhpb).eq.0.0d0)
1263      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1264         endif
1265         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1266      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1267         endif
1268 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1269 C        if (forcon(nhpb+1).gt.0.0d0) then
1270 C          nhpb=nhpb+1
1271 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1272       enddo
1273
1274       if (restr_type.eq.4) then
1275         write (iout,*) "The BFAC array"
1276         do i=nnt,nct
1277           write (iout,'(i5,f10.5)') i,bfac(i)
1278         enddo
1279         do i=nnt,nct
1280           if (itype(i).eq.ntyp1) cycle
1281           do j=nnt,i-1
1282             if (itype(j).eq.ntyp1) cycle
1283             if (itype(i).eq.10) then
1284               iiend=0
1285             else
1286               iiend=1
1287             endif
1288             if (itype(j).eq.10) then
1289               jjend=0
1290             else
1291               jjend=1
1292             endif
1293             kk=0
1294             do ii=0,iiend
1295             do jj=0,jjend
1296             nhpb=nhpb+1
1297             irestr_type(nhpb)=1
1298             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1299             irestr_type(nhpb)=1
1300             ibecarb(nhpb)=kk
1301             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1302             ihpb(nhpb)=i+nres*ii
1303             jhpb(nhpb)=j+nres*jj
1304             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1305             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1306      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1307      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1308      &       irestr_type(nhpb)
1309             kk=kk+1
1310           enddo
1311           enddo
1312           enddo
1313         enddo
1314       endif
1315
1316       if (restr_type.eq.5) then
1317         restr_on_coord=.true.
1318         do i=nnt,nct
1319           if (itype(i).eq.ntyp1) cycle
1320           bfac(i)=(scal_bfac/bfac(i))**2
1321         enddo
1322       endif
1323
1324       ENDDO ! next
1325
1326       fordepthmax=0.0d0
1327       if (normalize) then
1328         do i=nss+1,nhpb
1329           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1330      &      fordepthmax=fordepth(i)
1331         enddo
1332         do i=nss+1,nhpb
1333           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1334         enddo
1335       endif
1336       if (nhpb.gt.nss)  then
1337         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1338      &  "The following",nhpb-nss,
1339      &  " distance restraints have been imposed:",
1340      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
1341      &  "  score"," type"
1342         do i=nss+1,nhpb
1343           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1344      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1345      &  irestr_type(i)
1346         enddo
1347       endif
1348       call hpb_partition
1349       call flush(iout)
1350       return
1351    11 write (iout,*)"read_dist_restr: error reading reference structure"
1352       stop
1353       end
1354 c-------------------------------------------------------------------------------
1355       subroutine read_saxs_constr
1356       implicit real*8 (a-h,o-z)
1357       include 'DIMENSIONS'
1358 #ifdef MPI
1359       include 'mpif.h'
1360 #endif
1361       include 'COMMON.CONTROL'
1362       include 'COMMON.CHAIN'
1363       include 'COMMON.IOUNITS'
1364       include 'COMMON.SBRIDGE'
1365       include 'COMMON.SAXS'
1366       double precision cm(3)
1367 c      read(inp,*) nsaxs
1368       write (iout,*) "Calling read_saxs nsaxs",nsaxs
1369       call flush(iout)
1370       if (saxs_mode.eq.0) then
1371 c SAXS distance distribution
1372       do i=1,nsaxs
1373         read(inp,*) distsaxs(i),Psaxs(i)
1374       enddo
1375       Cnorm = 0.0d0
1376       do i=1,nsaxs
1377         Cnorm = Cnorm + Psaxs(i)
1378       enddo
1379       write (iout,*) "Cnorm",Cnorm
1380       do i=1,nsaxs
1381         Psaxs(i)=Psaxs(i)/Cnorm
1382       enddo
1383       write (iout,*) "Normalized distance distribution from SAXS"
1384       do i=1,nsaxs
1385         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1386       enddo
1387       Wsaxs0=0.0d0
1388       do i=1,nsaxs
1389         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1390       enddo
1391       write (iout,*) "Wsaxs0",Wsaxs0
1392       else
1393 c SAXS "spheres".
1394       do i=1,nsaxs
1395         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1396       enddo
1397       do j=1,3
1398         cm(j)=0.0d0
1399       enddo
1400       do i=1,nsaxs
1401         do j=1,3
1402           cm(j)=cm(j)+Csaxs(j,i)
1403         enddo
1404       enddo
1405       do j=1,3
1406         cm(j)=cm(j)/nsaxs
1407       enddo
1408       do i=1,nsaxs
1409         do j=1,3
1410           Csaxs(j,i)=Csaxs(j,i)-cm(j)
1411         enddo
1412       enddo
1413       write (iout,*) "SAXS sphere coordinates"
1414       do i=1,nsaxs
1415         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
1416       enddo
1417       endif
1418       return
1419       end