source/cluster/wham/src-M/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.INTERACT'
  17       include "COMMON.SPLITELE"
  18       include 'COMMON.SHIELD'
  19       include 'COMMON.FREE'
  20       include 'COMMON.LANGEVIN'
  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       call readi(controlcard,'NRES',nres,0)
  31       call readi(controlcard,'RESCALE',rescale_mode,2)
  32       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
  33       write (iout,*) "DISTCHAINMAX",distchainmax
  34 C Reading the dimensions of box in x,y,z coordinates
  35       call reada(controlcard,'BOXX',boxxsize,100.0d0)
  36       call reada(controlcard,'BOXY',boxysize,100.0d0)
  37       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
  38 c Cutoff range for interactions
  39       call reada(controlcard,"R_CUT",r_cut,15.0d0)
  40       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
  41       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
  42       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
  43       if (lipthick.gt.0.0d0) then
  44        bordliptop=(boxzsize+lipthick)/2.0
  45        bordlipbot=bordliptop-lipthick
  46 C      endif
  47       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
  48      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
  49       buflipbot=bordlipbot+lipbufthick
  50       bufliptop=bordliptop-lipbufthick
  51       if ((lipbufthick*2.0d0).gt.lipthick)
  52      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
  53       endif
  54       write(iout,*) "bordliptop=",bordliptop
  55       write(iout,*) "bordlipbot=",bordlipbot
  56       write(iout,*) "bufliptop=",bufliptop
  57       write(iout,*) "buflipbot=",buflipbot
  58 C Shielding mode
  59       call readi(controlcard,'SHIELD',shield_mode,0)
  60       write (iout,*) "SHIELD MODE",shield_mode
  61       if (shield_mode.gt.0) then
  62       pi=3.141592d0
  63 C VSolvSphere the volume of solving sphere
  64 C      print *,pi,"pi"
  65 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
  66 C there will be no distinction between proline peptide group and normal peptide
  67 C group in case of shielding parameters
  68       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
  69       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
  70       write (iout,*) VSolvSphere,VSolvSphere_div
  71 C long axis of side chain
  72       do i=1,ntyp
  73       long_r_sidechain(i)=vbldsc0(1,i)
  74       short_r_sidechain(i)=sigma0(i)
  75       enddo
  76       buff_shield=1.0d0
  77       endif
  78       call readi(controlcard,'PDBOUT',outpdb,0)
  79       call readi(controlcard,'MOL2OUT',outmol2,0)
  80       refstr=(index(controlcard,'REFSTR').gt.0)
  81       write (iout,*) "REFSTR",refstr
  82       pdbref=(index(controlcard,'PDBREF').gt.0)
  83       iscode=index(controlcard,'ONE_LETTER')
  84       tree=(index(controlcard,'MAKE_TREE').gt.0)
  85       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  86       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  87       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  88      & " CONSTR_DIST",constr_dist
  89       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
  90       write (iout,*) "with_theta_constr ",with_theta_constr
  91       call flush(iout)
  92       min_var=(index(controlcard,'MINVAR').gt.0)
  93       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
  94       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
  95       call readi(controlcard,'NCUT',ncut,0)
  96       call readi(controlcard,'NCLUST',nclust,5)
  97       call readi(controlcard,'SYM',symetr,1)
  98       write (iout,*) 'sym', symetr
  99       call readi(controlcard,'NSTART',nstart,0)
 100       call readi(controlcard,'NEND',nend,0)
 101       call reada(controlcard,'ECUT',ecut,10.0d0)
 102       call reada(controlcard,'PROB',prob_limit,0.99d0)
 103       write (iout,*) "Probability limit",prob_limit
 104       lgrp=(index(controlcard,'LGRP').gt.0)
 105       caonly=(index(controlcard,'CA_ONLY').gt.0)
 106       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
 107       if (ncut.gt.0)
 108      & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
 109       call readi(controlcard,'IOPT',iopt,2)
 110       lside = index(controlcard,"SIDE").gt.0
 111       efree = index(controlcard,"EFREE").gt.0
 112       call readi(controlcard,'NTEMP',nT,1)
 113       write (iout,*) "nT",nT
 114       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
 115       write (iout,*) "nT",nT
 116       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
 117       do i=1,nT
 118         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
 119       enddo
 120       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
 121       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
 122       lprint_int=index(controlcard,"PRINT_INT") .gt.0
 123       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
 124       write (iout,*) "with_homology_constr ",with_dihed_constr,
 125      & " CONSTR_HOMOLOGY",constr_homology
 126       print_homology_restraints=
 127      & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
 128       print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
 129       print_homology_models=
 130      & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
 131       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
 132       call readi(controlcard,'NSAXS',nsaxs,0)
 133       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 134       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 135       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 136       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 137      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 138       if (min_var) iopt=1
 139       return
 140       end
 141 c--------------------------------------------------------------------------
 142       subroutine molread
 143 C
 144 C Read molecular data.
 145 C
 146       implicit none
 147       include 'DIMENSIONS'
 148       include 'COMMON.IOUNITS'
 149       include 'COMMON.GEO'
 150       include 'COMMON.VAR'
 151       include 'COMMON.INTERACT'
 152       include 'COMMON.LOCAL'
 153       include 'COMMON.NAMES'
 154       include 'COMMON.CHAIN'
 155       include 'COMMON.FFIELD'
 156       include 'COMMON.SBRIDGE'
 157       include 'COMMON.HEADER'
 158       include 'COMMON.CONTROL'
 159       include 'COMMON.CONTACTS'
 160       include 'COMMON.TIME1'
 161       include 'COMMON.TORCNSTR'
 162       include 'COMMON.SHIELD'
 163 #ifdef MPL
 164       include 'COMMON.INFO'
 165 #endif
 166       character*4 sequence(maxres)
 167       character*800 weightcard
 168       integer rescode
 169       double precision x(maxvar)
 170       integer itype_pdb(maxres)
 171       logical seq_comp
 172       integer i,j,kkk,i1,i2,it1,it2
 173 C
 174 C Body
 175 C
 176 C Read weights of the subsequent energy terms.
 177       call card_concat(weightcard)
 178       write(iout,*) weightcard
 179       call reada(weightcard,'WSC',wsc,1.0d0)
 180       call reada(weightcard,'WLONG',wsc,wsc)
 181       call reada(weightcard,'WSCP',wscp,1.0d0)
 182       call reada(weightcard,'WELEC',welec,1.0D0)
 183       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 184       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 185       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 186       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 187       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 188       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 189       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 190       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 191       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 192       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 193       call reada(weightcard,'WBOND',wbond,1.0D0)
 194       call reada(weightcard,'WTOR',wtor,1.0D0)
 195       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 196       call reada(weightcard,'WANG',wang,1.0D0)
 197       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 198       call reada(weightcard,'WSAXS',wsaxs,0.0D0)
 199       call reada(weightcard,'SCAL14',scal14,0.4D0)
 200       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 201       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 202       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 203       if (index(weightcard,'SOFT').gt.0) ipot=6
 204       call reada(weightcard,"D0CM",d0cm,3.78d0)
 205       call reada(weightcard,"AKCM",akcm,15.1d0)
 206       call reada(weightcard,"AKTH",akth,11.0d0)
 207       call reada(weightcard,"AKCT",akct,12.0d0)
 208       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 209       call reada(weightcard,"V2SS",v2ss,7.61d0)
 210       call reada(weightcard,"V3SS",v3ss,13.7d0)
 211       call reada(weightcard,"EBR",ebr,-5.50D0)
 212       call reada(weightcard,'WSHIELD',wshield,1.0d0)
 213       write(iout,*) 'WSHIELD',wshield
 214        call reada(weightcard,'WLT',wliptran,0.0D0)
 215       call reada(weightcard,"ATRISS",atriss,0.301D0)
 216       call reada(weightcard,"BTRISS",btriss,0.021D0)
 217       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 218       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 219       write (iout,*) "ATRISS=", atriss
 220       write (iout,*) "BTRISS=", btriss
 221       write (iout,*) "CTRISS=", ctriss
 222       write (iout,*) "DTRISS=", dtriss
 223       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 224       do i=1,maxres
 225         dyn_ss_mask(i)=.false.
 226       enddo
 227       do i=1,maxres-1
 228         do j=i+1,maxres
 229           dyn_ssbond_ij(i,j)=1.0d300
 230         enddo
 231       enddo
 232       call reada(weightcard,"HT",Ht,0.0D0)
 233       if (dyn_ss) then
 234         ss_depth=ebr/wsc-0.25*eps(1,1)
 235         Ht=Ht/wsc-0.25*eps(1,1)
 236         akcm=akcm*wstrain/wsc
 237         akth=akth*wstrain/wsc
 238         akct=akct*wstrain/wsc
 239         v1ss=v1ss*wstrain/wsc
 240         v2ss=v2ss*wstrain/wsc
 241         v3ss=v3ss*wstrain/wsc
 242       else
 243         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 244       endif
 245       write (iout,'(/a)') "Disulfide bridge parameters:"
 246       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 247       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 248       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 249       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 250      & ' v3ss:',v3ss
 251
 252 C 12/1/95 Added weight for the multi-body term WCORR
 253       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 254       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 255       weights(1)=wsc
 256       weights(2)=wscp
 257       weights(3)=welec
 258       weights(4)=wcorr
 259       weights(5)=wcorr5
 260       weights(6)=wcorr6
 261       weights(7)=wel_loc
 262       weights(8)=wturn3
 263       weights(9)=wturn4
 264       weights(10)=wturn6
 265       weights(11)=wang
 266       weights(12)=wscloc
 267       weights(13)=wtor
 268       weights(14)=wtor_d
 269       weights(15)=wstrain
 270       weights(16)=wvdwpp
 271       weights(17)=wbond
 272       weights(18)=scal14
 273       weights(19)=wsccor
 274       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 275      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 276      &  wturn4,wturn6,wsccor
 277    10 format (/'Energy-term weights (unscaled):'//
 278      & 'WSCC=   ',f10.6,' (SC-SC)'/
 279      & 'WSCP=   ',f10.6,' (SC-p)'/
 280      & 'WELEC=  ',f10.6,' (p-p electr)'/
 281      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 282      & 'WBOND=  ',f10.6,' (stretching)'/
 283      & 'WANG=   ',f10.6,' (bending)'/
 284      & 'WSCLOC= ',f10.6,' (SC local)'/
 285      & 'WTOR=   ',f10.6,' (torsional)'/
 286      & 'WTORD=  ',f10.6,' (double torsional)'/
 287      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 288      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 289      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 290      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 291      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 292      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 293      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 294      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 295      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 296
 297       if (wcorr4.gt.0.0d0) then
 298         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 299      &   'between contact pairs of peptide groups'
 300         write (iout,'(2(a,f5.3/))')
 301      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 302      &  'Range of quenching the correlation terms:',2*delt_corr
 303       else if (wcorr.gt.0.0d0) then
 304         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 305      &   'between contact pairs of peptide groups'
 306       endif
 307       write (iout,'(a,f8.3)')
 308      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 309       write (iout,'(a,f8.3)')
 310      &  'General scaling factor of SC-p interactions:',scalscp
 311       r0_corr=cutoff_corr-delt_corr
 312       do i=1,20
 313         aad(i,1)=scalscp*aad(i,1)
 314         aad(i,2)=scalscp*aad(i,2)
 315         bad(i,1)=scalscp*bad(i,1)
 316         bad(i,2)=scalscp*bad(i,2)
 317       enddo
 318
 319       call flush(iout)
 320       print *,'indpdb=',indpdb,' pdbref=',pdbref
 321
 322 C Read sequence if not taken from the pdb file.
 323       if (iscode.gt.0) then
 324         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 325       else
 326         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 327       endif
 328 C Convert sequence to numeric code
 329       do i=1,nres
 330         itype(i)=rescode(i,sequence(i),iscode)
 331       enddo
 332       print *,nres
 333       print '(20i4)',(itype(i),i=1,nres)
 334
 335       do i=1,nres
 336 #ifdef PROCOR
 337         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 338 #else
 339         if (itype(i).eq.ntyp1) then
 340 #endif
 341           itel(i)=0
 342 #ifdef PROCOR
 343         else if (iabs(itype(i+1)).ne.20) then
 344 #else
 345         else if (iabs(itype(i)).ne.20) then
 346 #endif
 347           itel(i)=1
 348         else
 349           itel(i)=2
 350         endif
 351       enddo
 352       write (iout,*) "ITEL"
 353       do i=1,nres-1
 354         write (iout,*) i,itype(i),itel(i)
 355       enddo
 356
 357       print *,'Call Read_Bridge.'
 358       call read_bridge
 359 C this fragment reads diheadral constrains
 360       if (with_dihed_constr) then
 361
 362       read (inp,*) ndih_constr
 363       if (ndih_constr.gt.0) then
 364 C        read (inp,*) ftors
 365 C        write (iout,*) 'FTORS',ftors
 366 C ftors is the force constant for torsional quartic constrains
 367         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 368      &                i=1,ndih_constr)
 369         write (iout,*)
 370      &   'There are',ndih_constr,' constraints on phi angles.'
 371         do i=1,ndih_constr
 372           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 373      &  ftors(i)
 374         enddo
 375         do i=1,ndih_constr
 376           phi0(i)=deg2rad*phi0(i)
 377           drange(i)=deg2rad*drange(i)
 378         enddo
 379       endif ! endif ndif_constr.gt.0
 380       endif ! with_dihed_constr
 381       if (with_theta_constr) then
 382 C with_theta_constr is keyword allowing for occurance of theta constrains
 383       read (inp,*) ntheta_constr
 384 C ntheta_constr is the number of theta constrains
 385       if (ntheta_constr.gt.0) then
 386 C        read (inp,*) ftors
 387         read (inp,*) (itheta_constr(i),theta_constr0(i),
 388      &  theta_drange(i),for_thet_constr(i),
 389      &  i=1,ntheta_constr)
 390 C the above code reads from 1 to ntheta_constr
 391 C itheta_constr(i) residue i for which is theta_constr
 392 C theta_constr0 the global minimum value
 393 C theta_drange is range for which there is no energy penalty
 394 C for_thet_constr is the force constant for quartic energy penalty
 395 C E=k*x**4
 396 C        if(me.eq.king.or..not.out1file)then
 397          write (iout,*)
 398      &   'There are',ntheta_constr,' constraints on phi angles.'
 399          do i=1,ntheta_constr
 400           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 401      &    theta_drange(i),
 402      &    for_thet_constr(i)
 403          enddo
 404 C        endif
 405         do i=1,ntheta_constr
 406           theta_constr0(i)=deg2rad*theta_constr0(i)
 407           theta_drange(i)=deg2rad*theta_drange(i)
 408         enddo
 409 C        if(me.eq.king.or..not.out1file)
 410 C     &   write (iout,*) 'FTORS',ftors
 411 C        do i=1,ntheta_constr
 412 C          ii = itheta_constr(i)
 413 C          thetabound(1,ii) = phi0(i)-drange(i)
 414 C          thetabound(2,ii) = phi0(i)+drange(i)
 415 C        enddo
 416       endif ! ntheta_constr.gt.0
 417       endif! with_theta_constr
 418
 419       nnt=1
 420       nct=nres
 421       print *,'NNT=',NNT,' NCT=',NCT
 422       if (itype(1).eq.ntyp1) nnt=2
 423       if (itype(nres).eq.ntyp1) nct=nct-1
 424       if (nstart.lt.nnt) nstart=nnt
 425       if (nend.gt.nct .or. nend.eq.0) nend=nct
 426       write (iout,*) "nstart",nstart," nend",nend
 427       write (iout,*) "calling read_saxs_consrtr",nsaxs
 428       if (nsaxs.gt.0) call read_saxs_constr
 429       nres0=nres
 430       if (constr_homology.gt.0) then
 431         call read_constr_homology(print_homology_restraints)
 432       endif
 433
 434 c      if (pdbref) then
 435 c        read(inp,'(a)') pdbfile
 436 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 437 c        open(ipdbin,file=pdbfile,status='old',err=33)
 438 c        goto 34
 439 c  33    write (iout,'(a)') 'Error opening PDB file.'
 440 c        stop
 441 c  34    continue
 442 c        print *,'Begin reading pdb data'
 443 c        call readpdb
 444 c        print *,'Finished reading pdb data'
 445 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 446 c        do i=1,nres
 447 c          itype_pdb(i)=itype(i)
 448 c        enddo
 449 c        close (ipdbin)
 450 c        write (iout,'(a,i3)') 'nsup=',nsup
 451 c        nstart_seq=nnt
 452 c        if (nsup.le.(nct-nnt+1)) then
 453 c          do i=0,nct-nnt+1-nsup
 454 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 455 c              nstart_seq=nnt+i
 456 c              goto 111
 457 c            endif
 458 c          enddo
 459 c          write (iout,'(a)')
 460 c     &            'Error - sequences to be superposed do not match.'
 461 c          stop
 462 c        else
 463 c          do i=0,nsup-(nct-nnt+1)
 464 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 465 c     &      then
 466 c              nstart_sup=nstart_sup+i
 467 c              nsup=nct-nnt+1
 468 c              goto 111
 469 c            endif
 470 c          enddo
 471 c          write (iout,'(a)')
 472 c     &            'Error - sequences to be superposed do not match.'
 473 c        endif
 474 c  111   continue
 475 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 476 c     &                 ' nstart_seq=',nstart_seq
 477 c      endif
 478       call init_int_table
 479       call setup_var
 480       write (iout,*) "molread: REFSTR",refstr
 481       if (refstr) then
 482         if (.not.pdbref) then
 483           call read_angles(inp,*38)
 484           goto 39
 485    38     write (iout,'(a)') 'Error reading reference structure.'
 486 #ifdef MPL
 487           call mp_stopall(Error_Msg)
 488 #else
 489           stop 'Error reading reference structure'
 490 #endif
 491    39     call chainbuild
 492           nstart_sup=nnt
 493           nstart_seq=nnt
 494           nsup=nct-nnt+1
 495           kkk=1
 496           do i=1,2*nres
 497             do j=1,3
 498               cref(j,i,kkk)=c(j,i)
 499             enddo
 500           enddo
 501         endif
 502         call contact(.true.,ncont_ref,icont_ref)
 503       endif
 504        if (ns.gt.0) then
 505 C        write (iout,'(/a,i3,a)')
 506 C     &  'The chain contains',ns,' disulfide-bridging cysteines.'
 507         write (iout,'(20i4)') (iss(i),i=1,ns)
 508        if (dyn_ss) then
 509           write(iout,*)"Running with dynamic disulfide-bond formation"
 510        else
 511         write (iout,'(/a/)') 'Pre-formed links are:'
 512         do i=1,nss
 513           i1=ihpb(i)-nres
 514           i2=jhpb(i)-nres
 515           it1=itype(i1)
 516           it2=itype(i2)
 517           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 518      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 519      &    ebr,forcon(i)
 520         enddo
 521         write (iout,'(a)')
 522        endif
 523       endif
 524       if (ns.gt.0.and.dyn_ss) then
 525           do i=nss+1,nhpb
 526             ihpb(i-nss)=ihpb(i)
 527             jhpb(i-nss)=jhpb(i)
 528             forcon(i-nss)=forcon(i)
 529             dhpb(i-nss)=dhpb(i)
 530           enddo
 531           nhpb=nhpb-nss
 532           nss=0
 533           call hpb_partition
 534           do i=1,ns
 535             dyn_ss_mask(iss(i))=.true.
 536           enddo
 537       endif
 538 c Read distance restraints
 539       if (constr_dist.gt.0) then
 540         call read_dist_constr
 541         call hpb_partition
 542       endif
 543       return
 544       end
 545 c-----------------------------------------------------------------------------
 546       logical function seq_comp(itypea,itypeb,length)
 547       implicit none
 548       integer length,itypea(length),itypeb(length)
 549       integer i
 550       do i=1,length
 551         if (itypea(i).ne.itypeb(i)) then
 552           seq_comp=.false.
 553           return
 554         endif
 555       enddo
 556       seq_comp=.true.
 557       return
 558       end
 559 c-----------------------------------------------------------------------------
 560       subroutine read_bridge
 561 C Read information about disulfide bridges.
 562       implicit none
 563       include 'DIMENSIONS'
 564       include 'COMMON.IOUNITS'
 565       include 'COMMON.GEO'
 566       include 'COMMON.VAR'
 567       include 'COMMON.INTERACT'
 568       include 'COMMON.LOCAL'
 569       include 'COMMON.NAMES'
 570       include 'COMMON.CHAIN'
 571       include 'COMMON.FFIELD'
 572       include 'COMMON.SBRIDGE'
 573       include 'COMMON.HEADER'
 574       include 'COMMON.CONTROL'
 575       include 'COMMON.TIME1'
 576 #ifdef MPL
 577       include 'COMMON.INFO'
 578 #endif
 579       integer i,j
 580 C Read bridging residues.
 581       read (inp,*) ns,(iss(i),i=1,ns)
 582       print *,'ns=',ns
 583 C Check whether the specified bridging residues are cystines.
 584       do i=1,ns
 585         if (itype(iss(i)).ne.1) then
 586           write (iout,'(2a,i3,a)')
 587      &   'Do you REALLY think that the residue ',
 588      &    restyp(itype(iss(i))),i,
 589      &   ' can form a disulfide bridge?!!!'
 590           write (*,'(2a,i3,a)')
 591      &   'Do you REALLY think that the residue ',
 592      &   restyp(itype(iss(i))),i,
 593      &   ' can form a disulfide bridge?!!!'
 594 #ifdef MPL
 595          call mp_stopall(error_msg)
 596 #else
 597          stop
 598 #endif
 599         endif
 600       enddo
 601 C Read preformed bridges.
 602       if (ns.gt.0) then
 603       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 604       if (nss.gt.0) then
 605         nhpb=nss
 606 C Check if the residues involved in bridges are in the specified list of
 607 C bridging residues.
 608         do i=1,nss
 609           do j=1,i-1
 610             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 611      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 612               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 613      &      ' contains residues present in other pairs.'
 614               write (*,'(a,i3,a)') 'Disulfide pair',i,
 615      &      ' contains residues present in other pairs.'
 616 #ifdef MPL
 617               call mp_stopall(error_msg)
 618 #else
 619               stop
 620 #endif
 621             endif
 622           enddo
 623           do j=1,ns
 624             if (ihpb(i).eq.iss(j)) goto 10
 625           enddo
 626           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 627    10     continue
 628           do j=1,ns
 629             if (jhpb(i).eq.iss(j)) goto 20
 630           enddo
 631           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 632    20     continue
 633 C          dhpb(i)=dbr
 634 C          forcon(i)=fbr
 635         enddo
 636         do i=1,nss
 637           ihpb(i)=ihpb(i)+nres
 638           jhpb(i)=jhpb(i)+nres
 639         enddo
 640       endif
 641       endif
 642       return
 643       end
 644 c----------------------------------------------------------------------------
 645       subroutine read_angles(kanal,*)
 646       implicit none
 647       include 'DIMENSIONS'
 648       include 'COMMON.GEO'
 649       include 'COMMON.VAR'
 650       include 'COMMON.CHAIN'
 651       include 'COMMON.IOUNITS'
 652       integer i,kanal
 653       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 654       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 655       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 656       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 657       do i=1,nres
 658         theta(i)=deg2rad*theta(i)
 659         phi(i)=deg2rad*phi(i)
 660         alph(i)=deg2rad*alph(i)
 661         omeg(i)=deg2rad*omeg(i)
 662       enddo
 663       return
 664    10 return1
 665       end
 666 c----------------------------------------------------------------------------
 667       subroutine reada(rekord,lancuch,wartosc,default)
 668       implicit none
 669       character*(*) rekord,lancuch
 670       double precision wartosc,default
 671       integer ilen,iread
 672       external ilen
 673       iread=index(rekord,lancuch)
 674       if (iread.eq.0) then
 675         wartosc=default
 676         return
 677       endif
 678       iread=iread+ilen(lancuch)+1
 679       read (rekord(iread:),*) wartosc
 680       return
 681       end
 682 c----------------------------------------------------------------------------
 683       subroutine multreada(rekord,lancuch,tablica,dim,default)
 684       implicit none
 685       integer dim,i
 686       double precision tablica(dim),default
 687       character*(*) rekord,lancuch
 688       integer ilen,iread
 689       external ilen
 690       do i=1,dim
 691         tablica(i)=default
 692       enddo
 693       iread=index(rekord,lancuch)
 694       if (iread.eq.0) return
 695       iread=iread+ilen(lancuch)+1
 696       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 697    10 return
 698       end
 699 c----------------------------------------------------------------------------
 700       subroutine readi(rekord,lancuch,wartosc,default)
 701       implicit none
 702       character*(*) rekord,lancuch
 703       integer wartosc,default
 704       integer ilen,iread
 705       external ilen
 706       iread=index(rekord,lancuch)
 707       if (iread.eq.0) then
 708         wartosc=default
 709         return
 710       endif
 711       iread=iread+ilen(lancuch)+1
 712       read (rekord(iread:),*) wartosc
 713       return
 714       end
 715 C----------------------------------------------------------------------
 716       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 717       implicit none
 718       integer dim,i
 719       integer tablica(dim),default
 720       character*(*) rekord,lancuch
 721       character*80 aux
 722       integer ilen,iread
 723       external ilen
 724       do i=1,dim
 725         tablica(i)=default
 726       enddo
 727       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 728       if (iread.eq.0) return
 729       iread=iread+ilen(lancuch)+1
 730       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 731    10 return
 732       end
 733
 734 c----------------------------------------------------------------------------
 735       subroutine card_concat(card)
 736       include 'DIMENSIONS'
 737       include 'COMMON.IOUNITS'
 738       character*(*) card
 739       character*80 karta,ucase
 740       external ilen
 741       read (inp,'(a)') karta
 742       karta=ucase(karta)
 743       card=' '
 744       do while (karta(80:80).eq.'&')
 745         card=card(:ilen(card)+1)//karta(:79)
 746         read (inp,'(a)') karta
 747         karta=ucase(karta)
 748       enddo
 749       card=card(:ilen(card)+1)//karta
 750       return
 751       end
 752 c----------------------------------------------------------------------------
 753       subroutine openunits
 754       implicit none
 755       include 'DIMENSIONS'
 756 #ifdef MPI
 757       include "mpif.h"
 758       character*3 liczba
 759       include "COMMON.MPI"
 760 #endif
 761       include 'COMMON.IOUNITS'
 762       include 'COMMON.CONTROL'
 763       integer lenpre,lenpot,ilen
 764       external ilen
 765       character*16 cformat,cprint
 766       character*16 ucase
 767       integer lenint,lenout
 768       call getenv('INPUT',prefix)
 769       call getenv('OUTPUT',prefout)
 770       call getenv('INTIN',prefintin)
 771       call getenv('COORD',cformat)
 772       call getenv('PRINTCOOR',cprint)
 773       call getenv('SCRATCHDIR',scratchdir)
 774       from_bx=.true.
 775       from_cx=.false.
 776       if (index(ucase(cformat),'CX').gt.0) then
 777         from_cx=.true.
 778         from_bx=.false.
 779       endif
 780       from_cart=.true.
 781       lenpre=ilen(prefix)
 782       lenout=ilen(prefout)
 783       lenint=ilen(prefintin)
 784 C Get the names and open the input files
 785       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 786 #ifdef MPI
 787       write (liczba,'(bz,i3.3)') me
 788       outname=prefout(:lenout)//'_clust.out_'//liczba
 789 #else
 790       outname=prefout(:lenout)//'_clust.out'
 791 #endif
 792       if (from_bx) then
 793         intinname=prefintin(:lenint)//'.bx'
 794       else if (from_cx) then
 795         intinname=prefintin(:lenint)//'.cx'
 796       else
 797         intinname=prefintin(:lenint)//'.int'
 798       endif
 799       rmsname=prefintin(:lenint)//'.rms'
 800       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 801      &       status='unknown')
 802       open (jrms,file=rmsname,status='unknown')
 803       open(iout,file=outname,status='unknown')
 804 C Get parameter filenames and open the parameter files.
 805       call getenv('BONDPAR',bondname)
 806       open (ibond,file=bondname,status='old')
 807       call getenv('THETPAR',thetname)
 808       open (ithep,file=thetname,status='old')
 809       call getenv('ROTPAR',rotname)
 810       open (irotam,file=rotname,status='old')
 811       call getenv('TORPAR',torname)
 812       open (itorp,file=torname,status='old')
 813       call getenv('TORDPAR',tordname)
 814       open (itordp,file=tordname,status='old')
 815       call getenv('FOURIER',fouriername)
 816       open (ifourier,file=fouriername,status='old')
 817       call getenv('ELEPAR',elename)
 818       open (ielep,file=elename,status='old')
 819       call getenv('SIDEPAR',sidename)
 820       open (isidep,file=sidename,status='old')
 821       call getenv('SIDEP',sidepname)
 822       open (isidep1,file=sidepname,status="old")
 823       call getenv('SCCORPAR',sccorname)
 824       open (isccor,file=sccorname,status="old")
 825       call getenv('LIPTRANPAR',liptranname)
 826       open (iliptranpar,file=liptranname,status='old')
 827 #ifndef OLDSCP
 828 C
 829 C 8/9/01 In the newest version SCp interaction constants are read from a file
 830 C Use -DOLDSCP to use hard-coded constants instead.
 831 C
 832       call getenv('SCPPAR',scpname)
 833       open (iscpp,file=scpname,status='old')
 834 #endif
 835       return
 836       end
 837       subroutine read_dist_constr
 838       implicit real*8 (a-h,o-z)
 839       include 'DIMENSIONS'
 840 #ifdef MPI
 841       include 'mpif.h'
 842 #endif
 843       include 'COMMON.CONTROL'
 844       include 'COMMON.CHAIN'
 845       include 'COMMON.IOUNITS'
 846       include 'COMMON.SBRIDGE'
 847       integer ifrag_(2,100),ipair_(2,100)
 848       double precision wfrag_(100),wpair_(100)
 849       character*500 controlcard
 850       logical lprn /.true./
 851       logical normalize,next
 852       integer restr_type
 853       double precision xlink(4,0:4) /
 854 c           a          b       c     sigma
 855      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
 856      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
 857      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
 858      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
 859      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
 860       write (iout,*) "Calling read_dist_constr"
 861 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 862 c      call flush(iout)
 863       next=.true.
 864
 865       DO WHILE (next)
 866
 867       call card_concat(controlcard)
 868       next = index(controlcard,"NEXT").gt.0
 869       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
 870       write (iout,*) "restr_type",restr_type
 871       call readi(controlcard,"NFRAG",nfrag_,0)
 872       call readi(controlcard,"NFRAG",nfrag_,0)
 873       call readi(controlcard,"NPAIR",npair_,0)
 874       call readi(controlcard,"NDIST",ndist_,0)
 875       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 876       if (restr_type.eq.10)
 877      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
 878       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 879       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 880       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 881       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 882       normalize = index(controlcard,"NORMALIZE").gt.0
 883       write (iout,*) "WBOLTZD",wboltzd
 884 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 885 c      write (iout,*) "IFRAG"
 886 c      do i=1,nfrag_
 887 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 888 c      enddo
 889 c      write (iout,*) "IPAIR"
 890 c      do i=1,npair_
 891 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 892 c      enddo
 893       if (nfrag_.gt.0) then
 894         nres0=nres
 895         read(inp,'(a)') pdbfile
 896         write (iout,*)
 897      & "Distance restraints will be constructed from structure ",pdbfile
 898         open(ipdbin,file=pdbfile,status='old',err=11)
 899         call readpdb(.true.)
 900         nres=nres0
 901         close(ipdbin)
 902       endif
 903       do i=1,nfrag_
 904         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 905         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 906      &    ifrag_(2,i)=nstart_sup+nsup-1
 907 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 908 c        call flush(iout)
 909         if (wfrag_(i).eq.0.0d0) cycle
 910         do j=ifrag_(1,i),ifrag_(2,i)-1
 911           do k=j+1,ifrag_(2,i)
 912 c            write (iout,*) "j",j," k",k
 913             ddjk=dist(j,k)
 914             if (restr_type.eq.1) then
 915               nhpb=nhpb+1
 916               irestr_type(nhpb)=1
 917               ihpb(nhpb)=j
 918               jhpb(nhpb)=k
 919               dhpb(nhpb)=ddjk
 920               forcon(nhpb)=wfrag_(i)
 921             else if (constr_dist.eq.2) then
 922               if (ddjk.le.dist_cut) then
 923                 nhpb=nhpb+1
 924                 irestr_type(nhpb)=1
 925                 ihpb(nhpb)=j
 926                 jhpb(nhpb)=k
 927                 dhpb(nhpb)=ddjk
 928                 forcon(nhpb)=wfrag_(i)
 929               endif
 930             else if (restr_type.eq.3) then
 931               nhpb=nhpb+1
 932               irestr_type(nhpb)=1
 933               ihpb(nhpb)=j
 934               jhpb(nhpb)=k
 935               dhpb(nhpb)=ddjk
 936               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 937             endif
 938             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
 939      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 940           enddo
 941         enddo
 942       enddo
 943       do i=1,npair_
 944         if (wpair_(i).eq.0.0d0) cycle
 945         ii = ipair_(1,i)
 946         jj = ipair_(2,i)
 947         if (ii.gt.jj) then
 948           itemp=ii
 949           ii=jj
 950           jj=itemp
 951         endif
 952         do j=ifrag_(1,ii),ifrag_(2,ii)
 953           do k=ifrag_(1,jj),ifrag_(2,jj)
 954             if (restr_type.eq.1) then
 955               nhpb=nhpb+1
 956               irestr_type(nhpb)=1
 957               ihpb(nhpb)=j
 958               jhpb(nhpb)=k
 959               dhpb(nhpb)=ddjk
 960               forcon(nhpb)=wfrag_(i)
 961             else if (constr_dist.eq.2) then
 962               if (ddjk.le.dist_cut) then
 963                 nhpb=nhpb+1
 964                 irestr_type(nhpb)=1
 965                 ihpb(nhpb)=j
 966                 jhpb(nhpb)=k
 967                 dhpb(nhpb)=ddjk
 968                 forcon(nhpb)=wfrag_(i)
 969               endif
 970             else if (restr_type.eq.3) then
 971               nhpb=nhpb+1
 972               irestr_type(nhpb)=1
 973               ihpb(nhpb)=j
 974               jhpb(nhpb)=k
 975               dhpb(nhpb)=ddjk
 976               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 977             endif
 978             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
 979      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 980           enddo
 981         enddo
 982       enddo
 983
 984 c      print *,ndist_
 985       write (iout,*) "Distance restraints as read from input"
 986       do i=1,ndist_
 987         if (restr_type.eq.11) then
 988           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
 989      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 990 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 991           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
 992           nhpb=nhpb+1
 993           irestr_type(nhpb)=11
 994           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
 995      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 996      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 997           if (ibecarb(nhpb).gt.0) then
 998             ihpb(nhpb)=ihpb(nhpb)+nres
 999             jhpb(nhpb)=jhpb(nhpb)+nres
1000           endif
1001         else if (constr_dist.eq.10) then
1002 c Cross-lonk Markov-like potential
1003           call card_concat(controlcard)
1004           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1005           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1006           ibecarb(nhpb+1)=0
1007           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1008           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1009           if (index(controlcard,"ZL").gt.0) then
1010             link_type=1
1011           else if (index(controlcard,"ADH").gt.0) then
1012             link_type=2
1013           else if (index(controlcard,"PDH").gt.0) then
1014             link_type=3
1015           else if (index(controlcard,"DSS").gt.0) then
1016             link_type=4
1017           else
1018             link_type=0
1019           endif
1020           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1021      &       xlink(1,link_type))
1022           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1023      &       xlink(2,link_type))
1024           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1025      &       xlink(3,link_type))
1026           call reada(controlcard,"SIGMA",forcon(nhpb+1),
1027      &       xlink(4,link_type))
1028           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1029 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1030 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1031           if (forcon(nhpb+1).le.0.0d0 .or.
1032      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1033           nhpb=nhpb+1
1034           irestr_type(nhpb)=10
1035           if (ibecarb(nhpb).gt.0) then
1036             ihpb(nhpb)=ihpb(nhpb)+nres
1037             jhpb(nhpb)=jhpb(nhpb)+nres
1038           endif
1039           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1040      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1041      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1042      &     irestr_type(nhpb)
1043         else
1044 C        print *,"in else"
1045           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1046      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1047           if (forcon(nhpb+1).gt.0.0d0) then
1048           nhpb=nhpb+1
1049           if (dhpb1(nhpb).eq.0.0d0) then
1050             irestr_type(nhpb)=1
1051           else
1052             irestr_type(nhpb)=2
1053           endif
1054           if (ibecarb(nhpb).gt.0) then
1055             ihpb(nhpb)=ihpb(nhpb)+nres
1056             jhpb(nhpb)=jhpb(nhpb)+nres
1057           endif
1058           if (dhpb(nhpb).eq.0.0d0)
1059      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1060         endif
1061         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1062      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1063         endif
1064 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1065 C        if (forcon(nhpb+1).gt.0.0d0) then
1066 C          nhpb=nhpb+1
1067 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1068       enddo
1069
1070       ENDDO ! next
1071
1072       fordepthmax=0.0d0
1073       if (normalize) then
1074         do i=nss+1,nhpb
1075           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1076      &      fordepthmax=fordepth(i)
1077         enddo
1078         do i=nss+1,nhpb
1079           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1080         enddo
1081       endif
1082       if (nhpb.gt.nss)  then
1083         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1084      &  "The following",nhpb-nss,
1085      &  " distance restraints have been imposed:",
1086      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
1087      &  "  score"," type"
1088         do i=nss+1,nhpb
1089           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1090      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1091      &  irestr_type(i)
1092         enddo
1093       endif
1094       call hpb_partition
1095       call flush(iout)
1096       return
1097    11 write (iout,*)"read_dist_restr: error reading reference structure"
1098       stop
1099       end
1100 c-------------------------------------------------------------------------------
1101       subroutine read_saxs_constr
1102       implicit real*8 (a-h,o-z)
1103       include 'DIMENSIONS'
1104 #ifdef MPI
1105       include 'mpif.h'
1106 #endif
1107       include 'COMMON.CONTROL'
1108       include 'COMMON.CHAIN'
1109       include 'COMMON.IOUNITS'
1110       include 'COMMON.SBRIDGE'
1111       double precision cm(3)
1112 c      read(inp,*) nsaxs
1113       write (iout,*) "Calling read_saxs nsaxs",nsaxs
1114       call flush(iout)
1115       if (saxs_mode.eq.0) then
1116 c SAXS distance distribution
1117       do i=1,nsaxs
1118         read(inp,*) distsaxs(i),Psaxs(i)
1119       enddo
1120       Cnorm = 0.0d0
1121       do i=1,nsaxs
1122         Cnorm = Cnorm + Psaxs(i)
1123       enddo
1124       write (iout,*) "Cnorm",Cnorm
1125       do i=1,nsaxs
1126         Psaxs(i)=Psaxs(i)/Cnorm
1127       enddo
1128       write (iout,*) "Normalized distance distribution from SAXS"
1129       do i=1,nsaxs
1130         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1131       enddo
1132       Wsaxs0=0.0d0
1133       do i=1,nsaxs
1134         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1135       enddo
1136       write (iout,*) "Wsaxs0",Wsaxs0
1137       else
1138 c SAXS "spheres".
1139       do i=1,nsaxs
1140         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1141       enddo
1142       do j=1,3
1143         cm(j)=0.0d0
1144       enddo
1145       do i=1,nsaxs
1146         do j=1,3
1147           cm(j)=cm(j)+Csaxs(j,i)
1148         enddo
1149       enddo
1150       do j=1,3
1151         cm(j)=cm(j)/nsaxs
1152       enddo
1153       do i=1,nsaxs
1154         do j=1,3
1155           Csaxs(j,i)=Csaxs(j,i)-cm(j)
1156         enddo
1157       enddo
1158       write (iout,*) "SAXS sphere coordinates"
1159       do i=1,nsaxs
1160         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
1161       enddo
1162       endif
1163       return
1164       end
1165 c====-------------------------------------------------------------------
1166       subroutine read_constr_homology(lprn)
1167
1168       include 'DIMENSIONS'
1169 #ifdef MPI
1170       include 'mpif.h'
1171 #endif
1172       include 'COMMON.SETUP'
1173       include 'COMMON.CONTROL'
1174       include 'COMMON.CHAIN'
1175       include 'COMMON.IOUNITS'
1176       include 'COMMON.GEO'
1177       include 'COMMON.INTERACT'
1178       include 'COMMON.NAMES'
1179       include 'COMMON.HOMRESTR'
1180 c
1181 c For new homol impl
1182 c
1183       include 'COMMON.VAR'
1184 c     include 'include_unres/COMMON.VAR'
1185 c
1186
1187 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1188 c    &                 dist_cut
1189 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1190 c    &    sigma_odl_temp(maxres,maxres,max_template)
1191       character*2 kic2
1192       character*24 model_ki_dist, model_ki_angle
1193       character*500 controlcard
1194       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1195       integer idomain(max_template,maxres)
1196       integer ilen
1197       external ilen
1198       logical lprn
1199       logical unres_pdb,liiflag
1200 c
1201 c     FP - Nov. 2014 Temporary specifications for new vars
1202 c
1203       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
1204      &                 rescore3_tmp
1205       double precision, dimension (max_template,maxres) :: rescore
1206       double precision, dimension (max_template,maxres) :: rescore2
1207       double precision, dimension (max_template,maxres) :: rescore3
1208       character*24 tpl_k_rescore
1209 c -----------------------------------------------------------------
1210 c Reading multiple PDB ref structures and calculation of retraints
1211 c not using pre-computed ones stored in files model_ki_{dist,angle}
1212 c FP (Nov., 2014)
1213 c -----------------------------------------------------------------
1214 c
1215 c
1216 c Alternative: reading from input
1217 #ifdef DEBUG
1218       write (iout,*) "BEGIN READ HOMOLOGY INFO"
1219 #ifdef AIX
1220       call flush_(iout)
1221 #else
1222       call flush(iout)
1223 #endif
1224 #endif
1225       call card_concat(controlcard)
1226       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1227       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1228       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1229       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1230       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1231       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1232       dist1cut=(index(controlcard,'DIST1CUT').gt.0)
1233       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1234       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1235       if (homol_nset.gt.1)then
1236          call readi(controlcard,"ISET",iset,1)
1237          call card_concat(controlcard)
1238          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1239       else
1240         iset=1
1241         waga_homology(1)=1.0
1242       endif
1243 c
1244 #ifdef DEBUG
1245       write(iout,*) "read_constr_homology iset",iset
1246       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1247 #ifdef AIX
1248       call flush_(iout)
1249 #else
1250       call flush(iout)
1251 #endif
1252 #endif
1253 cd      write (iout,*) "nnt",nnt," nct",nct
1254 cd      call flush(iout)
1255
1256
1257       lim_odl=0
1258       lim_dih=0
1259 c
1260 c  New
1261 c
1262 c
1263 c  Reading HM global scores (prob not required)
1264 c
1265       do i = nnt,nct
1266         do k=1,constr_homology
1267          idomain(k,i)=0
1268         enddo
1269       enddo
1270 c     open (4,file="HMscore")
1271 c     do k=1,constr_homology
1272 c       read (4,*,end=521) hmscore_tmp
1273 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1274 c       write(*,*) "Model", k, ":", hmscore(k)
1275 c     enddo
1276 c521  continue
1277
1278       ii=0
1279       do i = nnt,nct-2
1280         do j=i+2,nct
1281         ii=ii+1
1282         ii_in_use(ii)=0
1283         enddo
1284       enddo
1285 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1286
1287       if (read_homol_frag) then
1288         call read_klapaucjusz
1289       else
1290       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1291       do k=1,constr_homology
1292
1293         read(inp,'(a)') pdbfile
1294 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1295 c  Next stament causes error upon compilation (?)
1296 c       if(me.eq.king.or. .not. out1file)
1297 c         write (iout,'(2a)') 'PDB data will be read from file ',
1298 c    &   pdbfile(:ilen(pdbfile))
1299          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1300      &  pdbfile(:ilen(pdbfile))
1301         open(ipdbin,file=pdbfile,status='old',err=33)
1302         goto 34
1303   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1304      &  pdbfile(:ilen(pdbfile))
1305         stop
1306   34    continue
1307 c        print *,'Begin reading pdb data'
1308 c
1309 c Files containing res sim or local scores (former containing sigmas)
1310 c
1311
1312         write(kic2,'(bz,i2.2)') k
1313
1314         tpl_k_rescore="template"//kic2//".sco"
1315
1316         unres_pdb=.false.
1317         call readpdb_template(k)
1318         do i=1,2*nres
1319           do j=1,3
1320             crefjlee(j,i)=c(j,i)
1321           enddo
1322         enddo
1323 #ifdef DEBUG
1324         do i=1,nres
1325           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1326      &      (crefjlee(j,i+nres),j=1,3)
1327         enddo
1328         write (iout,*) "READ HOMOLOGY INFO"
1329         write (iout,*) "read_constr_homology x: after reading pdb file"
1330         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1331         write (iout,*) "waga_dist",waga_dist
1332         write (iout,*) "waga_angle",waga_angle
1333         write (iout,*) "waga_theta",waga_theta
1334         write (iout,*) "waga_d",waga_d
1335         write (iout,*) "dist_cut",dist_cut
1336 #endif
1337 #ifdef AIX
1338       call flush_(iout)
1339 #else
1340       call flush(iout)
1341 #endif
1342
1343 c
1344 c     Distance restraints
1345 c
1346 c          ... --> odl(k,ii)
1347 C Copy the coordinates from reference coordinates (?)
1348         do i=1,2*nres
1349           do j=1,3
1350 c            c(j,i)=cref(j,i)
1351 c           write (iout,*) "c(",j,i,") =",c(j,i)
1352           enddo
1353         enddo
1354 c
1355 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1356 c
1357 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1358           open (ientin,file=tpl_k_rescore,status='old')
1359           if (nnt.gt.1) rescore(k,1)=0.0d0
1360           do irec=nnt,nct ! loop for reading res sim
1361             if (read2sigma) then
1362              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1363      &                                rescore3_tmp,idomain_tmp
1364              i_tmp=i_tmp+nnt-1
1365              idomain(k,i_tmp)=idomain_tmp
1366              rescore(k,i_tmp)=rescore_tmp
1367              rescore2(k,i_tmp)=rescore2_tmp
1368              rescore3(k,i_tmp)=rescore3_tmp
1369              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1370      &                      i_tmp,rescore2_tmp,rescore_tmp,
1371      &                                rescore3_tmp,idomain_tmp
1372             else
1373              idomain(k,irec)=1
1374              read (ientin,*,end=1401) rescore_tmp
1375
1376 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1377              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1378 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1379             endif
1380           enddo
1381  1401   continue
1382         close (ientin)
1383         if (waga_dist.ne.0.0d0) then
1384           ii=0
1385           do i = nnt,nct-2
1386             do j=i+2,nct
1387
1388               x12=c(1,i)-c(1,j)
1389               y12=c(2,i)-c(2,j)
1390               z12=c(3,i)-c(3,j)
1391               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1392 c              write (iout,*) k,i,j,distal,dist2_cut
1393            if (dist1cut .and. k.gt.1) then
1394               ii=ii+1
1395               if (l_homo(1,ii)) then
1396                 ii_in_use(ii)=1
1397                 l_homo(k,ii)=.true.
1398                 ires_homo(ii)=i
1399                 jres_homo(ii)=j
1400                 odl(k,ii)=distal
1401                 sigma_odl(k,ii)=sigma_odl(1,ii)
1402               else
1403                 l_homo(k,ii)=.false.
1404               endif
1405            else
1406             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1407      &            .and. distal.le.dist2_cut ) then
1408
1409               ii=ii+1
1410               ii_in_use(ii)=1
1411               l_homo(k,ii)=.true.
1412
1413 c             write (iout,*) "k",k
1414 c             write (iout,*) "i",i," j",j," constr_homology",
1415 c    &                       constr_homology
1416               ires_homo(ii)=i
1417               jres_homo(ii)=j
1418               odl(k,ii)=distal
1419               if (read2sigma) then
1420                 sigma_odl(k,ii)=0
1421                 do ik=i,j
1422                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1423                 enddo
1424                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1425                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1426      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1427               else
1428                 if (odl(k,ii).le.dist_cut) then
1429                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1430                 else
1431 #ifdef OLDSIGMA
1432                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1433      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1434 #else
1435                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1436      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1437 #endif
1438                 endif
1439               endif
1440               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1441             else
1442               ii=ii+1
1443               l_homo(k,ii)=.false.
1444             endif
1445            endif
1446             enddo
1447           enddo
1448         lim_odl=ii
1449         endif
1450 c
1451 c     Theta, dihedral and SC retraints
1452 c
1453         if (waga_angle.gt.0.0d0) then
1454 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1455 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1456 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1457 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1458 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1459 c    &                            sigma_dih(k,i+nnt-1)
1460 c         enddo
1461 c1402   continue
1462 c         close (ientin)
1463           do i = nnt+3,nct
1464             if (idomain(k,i).eq.0) then
1465                sigma_dih(k,i)=0.0
1466                cycle
1467             endif
1468             dih(k,i)=phiref(i) ! right?
1469 c           read (ientin,*) sigma_dih(k,i) ! original variant
1470 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1471 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1472 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1473 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1474 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1475
1476             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1477      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1478 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1479 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1480 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1481 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1482 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1483 c   Instead of res sim other local measure of b/b str reliability possible
1484             if (sigma_dih(k,i).ne.0)
1485      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1486 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1487           enddo
1488           lim_dih=nct-nnt-2
1489         endif
1490
1491         if (waga_theta.gt.0.0d0) then
1492 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1493 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1494 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1495 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1496 c    &                              sigma_theta(k,i+nnt-1)
1497 c         enddo
1498 c1403   continue
1499 c         close (ientin)
1500
1501           do i = nnt+2,nct ! right? without parallel.
1502 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1503 c         do i=ithet_start,ithet_end ! with FG parallel.
1504              if (idomain(k,i).eq.0) then
1505               sigma_theta(k,i)=0.0
1506               cycle
1507              endif
1508              thetatpl(k,i)=thetaref(i)
1509 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1510 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1511 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1512 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1513 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1514              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1515      &                        rescore(k,i-2))/3.0
1516 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1517              if (sigma_theta(k,i).ne.0)
1518      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1519
1520 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1521 c                             rescore(k,i-2) !  right expression ?
1522 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1523           enddo
1524         endif
1525
1526         if (waga_d.gt.0.0d0) then
1527 c       open (ientin,file=tpl_k_sigma_d,status='old')
1528 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1529 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1530 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1531 c    &                          sigma_d(k,i+nnt-1)
1532 c         enddo
1533 c1404   continue
1534
1535           do i = nnt,nct ! right? without parallel.
1536 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1537 c         do i=loc_start,loc_end ! with FG parallel.
1538                if (itype(i).eq.10) cycle
1539                if (idomain(k,i).eq.0 ) then
1540                   sigma_d(k,i)=0.0
1541                   cycle
1542                endif
1543                xxtpl(k,i)=xxref(i)
1544                yytpl(k,i)=yyref(i)
1545                zztpl(k,i)=zzref(i)
1546 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1547 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1548 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1549 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1550                sigma_d(k,i)=rescore3(k,i) !  right expression ?
1551                if (sigma_d(k,i).ne.0)
1552      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1553
1554 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1555 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1556 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1557           enddo
1558         endif
1559       enddo
1560 c
1561 c remove distance restraints not used in any model from the list
1562 c shift data in all arrays
1563 c
1564       if (waga_dist.ne.0.0d0) then
1565         ii=0
1566         liiflag=.true.
1567         do i=nnt,nct-2
1568          do j=i+2,nct
1569           ii=ii+1
1570           if (ii_in_use(ii).eq.0.and.liiflag) then
1571             liiflag=.false.
1572             iistart=ii
1573           endif
1574           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1575      &                   .not.liiflag.and.ii.eq.lim_odl) then
1576              if (ii.eq.lim_odl) then
1577               iishift=ii-iistart+1
1578              else
1579               iishift=ii-iistart
1580              endif
1581              liiflag=.true.
1582              do ki=iistart,lim_odl-iishift
1583               ires_homo(ki)=ires_homo(ki+iishift)
1584               jres_homo(ki)=jres_homo(ki+iishift)
1585               ii_in_use(ki)=ii_in_use(ki+iishift)
1586               do k=1,constr_homology
1587                odl(k,ki)=odl(k,ki+iishift)
1588                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1589                l_homo(k,ki)=l_homo(k,ki+iishift)
1590               enddo
1591              enddo
1592              ii=ii-iishift
1593              lim_odl=lim_odl-iishift
1594           endif
1595          enddo
1596         enddo
1597       endif
1598
1599       endif ! .not. klapaucjusz
1600
1601       if (constr_homology.gt.0) call homology_partition
1602       if (constr_homology.gt.0) call init_int_table
1603 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1604 cd     & "lim_xx=",lim_xx
1605 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1606 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1607 c
1608 c Print restraints
1609 c
1610       if (.not.lprn) return
1611 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1612       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1613        write (iout,*) "Distance restraints from templates"
1614        do ii=1,lim_odl
1615        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1616      &  ii,ires_homo(ii),jres_homo(ii),
1617      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1618      &  ki=1,constr_homology)
1619        enddo
1620        write (iout,*) "Dihedral angle restraints from templates"
1621        do i=nnt+3,nct
1622         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1623      &      (rad2deg*dih(ki,i),
1624      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1625        enddo
1626        write (iout,*) "Virtual-bond angle restraints from templates"
1627        do i=nnt+2,nct
1628         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1629      &      (rad2deg*thetatpl(ki,i),
1630      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1631        enddo
1632        write (iout,*) "SC restraints from templates"
1633        do i=nnt,nct
1634         write(iout,'(i5,100(4f8.2,4x))') i,
1635      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1636      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1637        enddo
1638       endif
1639 c -----------------------------------------------------------------
1640       return
1641       end
1642 c----------------------------------------------------------------------
1643       subroutine read_klapaucjusz
1644
1645       include 'DIMENSIONS'
1646 #ifdef MPI
1647       include 'mpif.h'
1648 #endif
1649       include 'COMMON.SETUP'
1650       include 'COMMON.CONTROL'
1651       include 'COMMON.CHAIN'
1652       include 'COMMON.IOUNITS'
1653       include 'COMMON.GEO'
1654       include 'COMMON.INTERACT'
1655       include 'COMMON.NAMES'
1656       include 'COMMON.HOMRESTR'
1657       character*256 fragfile
1658       integer ninclust(maxclust),inclust(max_template,maxclust),
1659      &  nresclust(maxclust),iresclust(maxres,maxclust)
1660
1661       character*2 kic2
1662       character*24 model_ki_dist, model_ki_angle
1663       character*500 controlcard
1664       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1665       integer idomain(max_template,maxres)
1666       logical lprn /.true./
1667       integer ilen
1668       external ilen
1669       logical unres_pdb,liiflag
1670 c
1671 c
1672       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1673       double precision, dimension (max_template,maxres) :: rescore
1674       double precision, dimension (max_template,maxres) :: rescore2
1675       character*24 tpl_k_rescore
1676
1677 c
1678 c For new homol impl
1679 c
1680       include 'COMMON.VAR'
1681 c
1682       double precision chomo(3,maxres2+2,max_template)
1683       call getenv("FRAGFILE",fragfile)
1684       open(ientin,file=fragfile,status="old",err=10)
1685       read(ientin,*) constr_homology,nclust
1686       l_homo = .false.
1687       sigma_theta=0.0
1688       sigma_d=0.0
1689       sigma_dih=0.0
1690 c Read pdb files
1691       do k=1,constr_homology
1692         read(ientin,'(a)') pdbfile
1693         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1694      &  pdbfile(:ilen(pdbfile))
1695         open(ipdbin,file=pdbfile,status='old',err=33)
1696         goto 34
1697   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1698      &  pdbfile(:ilen(pdbfile))
1699         stop
1700   34    continue
1701         unres_pdb=.false.
1702         call readpdb_template(k)
1703         do i=1,2*nres
1704           do j=1,3
1705             chomo(j,i,k)=c(j,i)
1706           enddo
1707         enddo
1708         do i=1,nres
1709           rescore(k,i)=0.2d0
1710           rescore2(k,i)=1.0d0
1711         enddo
1712       enddo
1713 c Read clusters
1714       do i=1,nclust
1715         read(ientin,*) ninclust(i),nresclust(i)
1716         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1717         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1718       enddo
1719 c
1720 c Loop over clusters
1721 c
1722       do l=1,nclust
1723         do ll = 1,ninclust(l)
1724
1725         k = inclust(ll,l)
1726         do i=1,nres
1727           idomain(k,i)=0
1728         enddo
1729         do i=1,nresclust(l)
1730           if (nnt.gt.1)  then
1731             idomain(k,iresclust(i,l)+1) = 1
1732           else
1733             idomain(k,iresclust(i,l)) = 1
1734           endif
1735         enddo
1736 c
1737 c     Distance restraints
1738 c
1739 c          ... --> odl(k,ii)
1740 C Copy the coordinates from reference coordinates (?)
1741         do i=1,2*nres
1742           do j=1,3
1743             c(j,i)=chomo(j,i,k)
1744 c           write (iout,*) "c(",j,i,") =",c(j,i)
1745           enddo
1746         enddo
1747         call int_from_cart1(.false.)
1748         call int_from_cart(.true.,.false.)
1749         call sc_loc_geom(.false.)
1750         do i=1,nres
1751           thetaref(i)=theta(i)
1752           phiref(i)=phi(i)
1753         enddo
1754         if (waga_dist.ne.0.0d0) then
1755           ii=0
1756           do i = nnt,nct-2
1757             do j=i+2,nct
1758
1759               x12=c(1,i)-c(1,j)
1760               y12=c(2,i)-c(2,j)
1761               z12=c(3,i)-c(3,j)
1762               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1763 c              write (iout,*) k,i,j,distal,dist2_cut
1764
1765             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1766      &            .and. distal.le.dist2_cut ) then
1767
1768               ii=ii+1
1769               ii_in_use(ii)=1
1770               l_homo(k,ii)=.true.
1771
1772 c             write (iout,*) "k",k
1773 c             write (iout,*) "i",i," j",j," constr_homology",
1774 c    &                       constr_homology
1775               ires_homo(ii)=i
1776               jres_homo(ii)=j
1777               odl(k,ii)=distal
1778               if (read2sigma) then
1779                 sigma_odl(k,ii)=0
1780                 do ik=i,j
1781                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1782                 enddo
1783                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1784                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1785      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1786               else
1787                 if (odl(k,ii).le.dist_cut) then
1788                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1789                 else
1790 #ifdef OLDSIGMA
1791                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1792      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1793 #else
1794                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1795      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1796 #endif
1797                 endif
1798               endif
1799               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1800             else
1801               ii=ii+1
1802 c              l_homo(k,ii)=.false.
1803             endif
1804             enddo
1805           enddo
1806         lim_odl=ii
1807         endif
1808 c
1809 c     Theta, dihedral and SC retraints
1810 c
1811         if (waga_angle.gt.0.0d0) then
1812           do i = nnt+3,nct
1813             if (idomain(k,i).eq.0) then
1814 c               sigma_dih(k,i)=0.0
1815                cycle
1816             endif
1817             dih(k,i)=phiref(i)
1818             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1819      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1820 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1821 c     &       " sigma_dihed",sigma_dih(k,i)
1822             if (sigma_dih(k,i).ne.0)
1823      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1824           enddo
1825           lim_dih=nct-nnt-2
1826         endif
1827
1828         if (waga_theta.gt.0.0d0) then
1829           do i = nnt+2,nct
1830              if (idomain(k,i).eq.0) then
1831 c              sigma_theta(k,i)=0.0
1832               cycle
1833              endif
1834              thetatpl(k,i)=thetaref(i)
1835              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1836      &                        rescore(k,i-2))/3.0
1837              if (sigma_theta(k,i).ne.0)
1838      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1839           enddo
1840         endif
1841
1842         if (waga_d.gt.0.0d0) then
1843           do i = nnt,nct
1844                if (itype(i).eq.10) cycle
1845                if (idomain(k,i).eq.0 ) then
1846 c                  sigma_d(k,i)=0.0
1847                   cycle
1848                endif
1849                xxtpl(k,i)=xxref(i)
1850                yytpl(k,i)=yyref(i)
1851                zztpl(k,i)=zzref(i)
1852                sigma_d(k,i)=rescore(k,i)
1853                if (sigma_d(k,i).ne.0)
1854      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1855                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1856           enddo
1857         endif
1858       enddo ! l
1859       enddo ! ll
1860 c
1861 c remove distance restraints not used in any model from the list
1862 c shift data in all arrays
1863 c
1864       if (waga_dist.ne.0.0d0) then
1865         ii=0
1866         liiflag=.true.
1867         do i=nnt,nct-2
1868          do j=i+2,nct
1869           ii=ii+1
1870           if (ii_in_use(ii).eq.0.and.liiflag) then
1871             liiflag=.false.
1872             iistart=ii
1873           endif
1874           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1875      &                   .not.liiflag.and.ii.eq.lim_odl) then
1876              if (ii.eq.lim_odl) then
1877               iishift=ii-iistart+1
1878              else
1879               iishift=ii-iistart
1880              endif
1881              liiflag=.true.
1882              do ki=iistart,lim_odl-iishift
1883               ires_homo(ki)=ires_homo(ki+iishift)
1884               jres_homo(ki)=jres_homo(ki+iishift)
1885               ii_in_use(ki)=ii_in_use(ki+iishift)
1886               do k=1,constr_homology
1887                odl(k,ki)=odl(k,ki+iishift)
1888                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1889                l_homo(k,ki)=l_homo(k,ki+iishift)
1890               enddo
1891              enddo
1892              ii=ii-iishift
1893              lim_odl=lim_odl-iishift
1894           endif
1895          enddo
1896         enddo
1897       endif
1898
1899       return
1900    10 stop "Error infragment file"
1901       end
1902 c----------------------------------------------------------------------