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       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1233       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1234       if (homol_nset.gt.1)then
1235          call readi(controlcard,"ISET",iset,1)
1236          call card_concat(controlcard)
1237          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1238       else
1239         iset=1
1240         waga_homology(1)=1.0
1241       endif
1242 c
1243 #ifdef DEBUG
1244       write(iout,*) "read_constr_homology iset",iset
1245       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1246 #ifdef AIX
1247       call flush_(iout)
1248 #else
1249       call flush(iout)
1250 #endif
1251 #endif
1252 cd      write (iout,*) "nnt",nnt," nct",nct
1253 cd      call flush(iout)
1254
1255
1256       lim_odl=0
1257       lim_dih=0
1258 c
1259 c  New
1260 c
1261 c
1262 c  Reading HM global scores (prob not required)
1263 c
1264       do i = nnt,nct
1265         do k=1,constr_homology
1266          idomain(k,i)=0
1267         enddo
1268       enddo
1269 c     open (4,file="HMscore")
1270 c     do k=1,constr_homology
1271 c       read (4,*,end=521) hmscore_tmp
1272 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1273 c       write(*,*) "Model", k, ":", hmscore(k)
1274 c     enddo
1275 c521  continue
1276
1277       ii=0
1278       do i = nnt,nct-2
1279         do j=i+2,nct
1280         ii=ii+1
1281         ii_in_use(ii)=0
1282         enddo
1283       enddo
1284 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1285
1286       if (read_homol_frag) then
1287         call read_klapaucjusz
1288       else
1289       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1290       do k=1,constr_homology
1291
1292         read(inp,'(a)') pdbfile
1293 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1294 c  Next stament causes error upon compilation (?)
1295 c       if(me.eq.king.or. .not. out1file)
1296 c         write (iout,'(2a)') 'PDB data will be read from file ',
1297 c    &   pdbfile(:ilen(pdbfile))
1298          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1299      &  pdbfile(:ilen(pdbfile))
1300         open(ipdbin,file=pdbfile,status='old',err=33)
1301         goto 34
1302   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1303      &  pdbfile(:ilen(pdbfile))
1304         stop
1305   34    continue
1306 c        print *,'Begin reading pdb data'
1307 c
1308 c Files containing res sim or local scores (former containing sigmas)
1309 c
1310
1311         write(kic2,'(bz,i2.2)') k
1312
1313         tpl_k_rescore="template"//kic2//".sco"
1314
1315         unres_pdb=.false.
1316         call readpdb_template(k)
1317         do i=1,2*nres
1318           do j=1,3
1319             crefjlee(j,i)=c(j,i)
1320           enddo
1321         enddo
1322 #ifdef DEBUG
1323         do i=1,nres
1324           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1325      &      (crefjlee(j,i+nres),j=1,3)
1326         enddo
1327         write (iout,*) "READ HOMOLOGY INFO"
1328         write (iout,*) "read_constr_homology x: after reading pdb file"
1329         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1330         write (iout,*) "waga_dist",waga_dist
1331         write (iout,*) "waga_angle",waga_angle
1332         write (iout,*) "waga_theta",waga_theta
1333         write (iout,*) "waga_d",waga_d
1334         write (iout,*) "dist_cut",dist_cut
1335 #endif
1336 #ifdef AIX
1337       call flush_(iout)
1338 #else
1339       call flush(iout)
1340 #endif
1341
1342 c
1343 c     Distance restraints
1344 c
1345 c          ... --> odl(k,ii)
1346 C Copy the coordinates from reference coordinates (?)
1347         do i=1,2*nres
1348           do j=1,3
1349 c            c(j,i)=cref(j,i)
1350 c           write (iout,*) "c(",j,i,") =",c(j,i)
1351           enddo
1352         enddo
1353 c
1354 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1355 c
1356 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1357           open (ientin,file=tpl_k_rescore,status='old')
1358           if (nnt.gt.1) rescore(k,1)=0.0d0
1359           do irec=nnt,nct ! loop for reading res sim
1360             if (read2sigma) then
1361              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1362      &                                rescore3_tmp,idomain_tmp
1363              i_tmp=i_tmp+nnt-1
1364              idomain(k,i_tmp)=idomain_tmp
1365              rescore(k,i_tmp)=rescore_tmp
1366              rescore2(k,i_tmp)=rescore2_tmp
1367              rescore3(k,i_tmp)=rescore3_tmp
1368              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1369      &                      i_tmp,rescore2_tmp,rescore_tmp,
1370      &                                rescore3_tmp,idomain_tmp
1371             else
1372              idomain(k,irec)=1
1373              read (ientin,*,end=1401) rescore_tmp
1374
1375 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1376              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1377 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1378             endif
1379           enddo
1380  1401   continue
1381         close (ientin)
1382         if (waga_dist.ne.0.0d0) then
1383           ii=0
1384           do i = nnt,nct-2
1385             do j=i+2,nct
1386
1387               x12=c(1,i)-c(1,j)
1388               y12=c(2,i)-c(2,j)
1389               z12=c(3,i)-c(3,j)
1390               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1391 c              write (iout,*) k,i,j,distal,dist2_cut
1392
1393             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1394      &            .and. distal.le.dist2_cut ) then
1395
1396               ii=ii+1
1397               ii_in_use(ii)=1
1398               l_homo(k,ii)=.true.
1399
1400 c             write (iout,*) "k",k
1401 c             write (iout,*) "i",i," j",j," constr_homology",
1402 c    &                       constr_homology
1403               ires_homo(ii)=i
1404               jres_homo(ii)=j
1405               odl(k,ii)=distal
1406               if (read2sigma) then
1407                 sigma_odl(k,ii)=0
1408                 do ik=i,j
1409                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1410                 enddo
1411                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1412                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1413      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1414               else
1415                 if (odl(k,ii).le.dist_cut) then
1416                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1417                 else
1418 #ifdef OLDSIGMA
1419                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1420      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1421 #else
1422                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1423      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1424 #endif
1425                 endif
1426               endif
1427               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1428             else
1429               ii=ii+1
1430               l_homo(k,ii)=.false.
1431             endif
1432             enddo
1433           enddo
1434         lim_odl=ii
1435         endif
1436 c
1437 c     Theta, dihedral and SC retraints
1438 c
1439         if (waga_angle.gt.0.0d0) then
1440 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1441 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1442 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1443 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1444 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1445 c    &                            sigma_dih(k,i+nnt-1)
1446 c         enddo
1447 c1402   continue
1448 c         close (ientin)
1449           do i = nnt+3,nct
1450             if (idomain(k,i).eq.0) then
1451                sigma_dih(k,i)=0.0
1452                cycle
1453             endif
1454             dih(k,i)=phiref(i) ! right?
1455 c           read (ientin,*) sigma_dih(k,i) ! original variant
1456 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1457 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1458 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1459 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1460 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1461
1462             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1463      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1464 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1465 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1466 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1467 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1468 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1469 c   Instead of res sim other local measure of b/b str reliability possible
1470             if (sigma_dih(k,i).ne.0)
1471      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1472 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1473           enddo
1474           lim_dih=nct-nnt-2
1475         endif
1476
1477         if (waga_theta.gt.0.0d0) then
1478 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1479 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1480 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1481 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1482 c    &                              sigma_theta(k,i+nnt-1)
1483 c         enddo
1484 c1403   continue
1485 c         close (ientin)
1486
1487           do i = nnt+2,nct ! right? without parallel.
1488 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1489 c         do i=ithet_start,ithet_end ! with FG parallel.
1490              if (idomain(k,i).eq.0) then
1491               sigma_theta(k,i)=0.0
1492               cycle
1493              endif
1494              thetatpl(k,i)=thetaref(i)
1495 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1496 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1497 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1498 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1499 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1500              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1501      &                        rescore(k,i-2))/3.0
1502 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1503              if (sigma_theta(k,i).ne.0)
1504      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1505
1506 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1507 c                             rescore(k,i-2) !  right expression ?
1508 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1509           enddo
1510         endif
1511
1512         if (waga_d.gt.0.0d0) then
1513 c       open (ientin,file=tpl_k_sigma_d,status='old')
1514 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1515 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1516 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1517 c    &                          sigma_d(k,i+nnt-1)
1518 c         enddo
1519 c1404   continue
1520
1521           do i = nnt,nct ! right? without parallel.
1522 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1523 c         do i=loc_start,loc_end ! with FG parallel.
1524                if (itype(i).eq.10) cycle
1525                if (idomain(k,i).eq.0 ) then
1526                   sigma_d(k,i)=0.0
1527                   cycle
1528                endif
1529                xxtpl(k,i)=xxref(i)
1530                yytpl(k,i)=yyref(i)
1531                zztpl(k,i)=zzref(i)
1532 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1533 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1534 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1535 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1536                sigma_d(k,i)=rescore3(k,i) !  right expression ?
1537                if (sigma_d(k,i).ne.0)
1538      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1539
1540 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1541 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1542 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1543           enddo
1544         endif
1545       enddo
1546 c
1547 c remove distance restraints not used in any model from the list
1548 c shift data in all arrays
1549 c
1550       if (waga_dist.ne.0.0d0) then
1551         ii=0
1552         liiflag=.true.
1553         do i=nnt,nct-2
1554          do j=i+2,nct
1555           ii=ii+1
1556           if (ii_in_use(ii).eq.0.and.liiflag) then
1557             liiflag=.false.
1558             iistart=ii
1559           endif
1560           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1561      &                   .not.liiflag.and.ii.eq.lim_odl) then
1562              if (ii.eq.lim_odl) then
1563               iishift=ii-iistart+1
1564              else
1565               iishift=ii-iistart
1566              endif
1567              liiflag=.true.
1568              do ki=iistart,lim_odl-iishift
1569               ires_homo(ki)=ires_homo(ki+iishift)
1570               jres_homo(ki)=jres_homo(ki+iishift)
1571               ii_in_use(ki)=ii_in_use(ki+iishift)
1572               do k=1,constr_homology
1573                odl(k,ki)=odl(k,ki+iishift)
1574                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1575                l_homo(k,ki)=l_homo(k,ki+iishift)
1576               enddo
1577              enddo
1578              ii=ii-iishift
1579              lim_odl=lim_odl-iishift
1580           endif
1581          enddo
1582         enddo
1583       endif
1584
1585       endif ! .not. klapaucjusz
1586
1587       if (constr_homology.gt.0) call homology_partition
1588       if (constr_homology.gt.0) call init_int_table
1589 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1590 cd     & "lim_xx=",lim_xx
1591 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1592 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1593 c
1594 c Print restraints
1595 c
1596       if (.not.lprn) return
1597 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1598       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1599        write (iout,*) "Distance restraints from templates"
1600        do ii=1,lim_odl
1601        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1602      &  ii,ires_homo(ii),jres_homo(ii),
1603      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1604      &  ki=1,constr_homology)
1605        enddo
1606        write (iout,*) "Dihedral angle restraints from templates"
1607        do i=nnt+3,nct
1608         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1609      &      (rad2deg*dih(ki,i),
1610      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1611        enddo
1612        write (iout,*) "Virtual-bond angle restraints from templates"
1613        do i=nnt+2,nct
1614         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1615      &      (rad2deg*thetatpl(ki,i),
1616      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1617        enddo
1618        write (iout,*) "SC restraints from templates"
1619        do i=nnt,nct
1620         write(iout,'(i5,100(4f8.2,4x))') i,
1621      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1622      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1623        enddo
1624       endif
1625 c -----------------------------------------------------------------
1626       return
1627       end
1628 c----------------------------------------------------------------------
1629       subroutine read_klapaucjusz
1630
1631       include 'DIMENSIONS'
1632 #ifdef MPI
1633       include 'mpif.h'
1634 #endif
1635       include 'COMMON.SETUP'
1636       include 'COMMON.CONTROL'
1637       include 'COMMON.CHAIN'
1638       include 'COMMON.IOUNITS'
1639       include 'COMMON.GEO'
1640       include 'COMMON.INTERACT'
1641       include 'COMMON.NAMES'
1642       include 'COMMON.HOMRESTR'
1643       character*256 fragfile
1644       integer ninclust(maxclust),inclust(max_template,maxclust),
1645      &  nresclust(maxclust),iresclust(maxres,maxclust)
1646
1647       character*2 kic2
1648       character*24 model_ki_dist, model_ki_angle
1649       character*500 controlcard
1650       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1651       integer idomain(max_template,maxres)
1652       logical lprn /.true./
1653       integer ilen
1654       external ilen
1655       logical unres_pdb,liiflag
1656 c
1657 c
1658       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1659       double precision, dimension (max_template,maxres) :: rescore
1660       double precision, dimension (max_template,maxres) :: rescore2
1661       character*24 tpl_k_rescore
1662
1663 c
1664 c For new homol impl
1665 c
1666       include 'COMMON.VAR'
1667 c
1668       double precision chomo(3,maxres2+2,max_template)
1669       call getenv("FRAGFILE",fragfile)
1670       open(ientin,file=fragfile,status="old",err=10)
1671       read(ientin,*) constr_homology,nclust
1672       l_homo = .false.
1673       sigma_theta=0.0
1674       sigma_d=0.0
1675       sigma_dih=0.0
1676 c Read pdb files
1677       do k=1,constr_homology
1678         read(ientin,'(a)') pdbfile
1679         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1680      &  pdbfile(:ilen(pdbfile))
1681         open(ipdbin,file=pdbfile,status='old',err=33)
1682         goto 34
1683   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1684      &  pdbfile(:ilen(pdbfile))
1685         stop
1686   34    continue
1687         unres_pdb=.false.
1688         call readpdb
1689         do i=1,2*nres
1690           do j=1,3
1691             chomo(j,i,k)=c(j,i)
1692           enddo
1693         enddo
1694         do i=1,nres
1695           rescore(k,i)=0.2d0
1696           rescore2(k,i)=1.0d0
1697         enddo
1698       enddo
1699 c Read clusters
1700       do i=1,nclust
1701         read(ientin,*) ninclust(i),nresclust(i)
1702         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1703         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1704       enddo
1705 c
1706 c Loop over clusters
1707 c
1708       do l=1,nclust
1709         do ll = 1,ninclust(l)
1710
1711         k = inclust(ll,l)
1712         do i=1,nres
1713           idomain(k,i)=0
1714         enddo
1715         do i=1,nresclust(l)
1716           if (nnt.gt.1)  then
1717             idomain(k,iresclust(i,l)+1) = 1
1718           else
1719             idomain(k,iresclust(i,l)) = 1
1720           endif
1721         enddo
1722 c
1723 c     Distance restraints
1724 c
1725 c          ... --> odl(k,ii)
1726 C Copy the coordinates from reference coordinates (?)
1727         do i=1,2*nres
1728           do j=1,3
1729             c(j,i)=chomo(j,i,k)
1730 c           write (iout,*) "c(",j,i,") =",c(j,i)
1731           enddo
1732         enddo
1733         call int_from_cart(.true.,.false.)
1734         call sc_loc_geom(.false.)
1735         do i=1,nres
1736           thetaref(i)=theta(i)
1737           phiref(i)=phi(i)
1738         enddo
1739         if (waga_dist.ne.0.0d0) then
1740           ii=0
1741           do i = nnt,nct-2
1742             do j=i+2,nct
1743
1744               x12=c(1,i)-c(1,j)
1745               y12=c(2,i)-c(2,j)
1746               z12=c(3,i)-c(3,j)
1747               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1748 c              write (iout,*) k,i,j,distal,dist2_cut
1749
1750             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1751      &            .and. distal.le.dist2_cut ) then
1752
1753               ii=ii+1
1754               ii_in_use(ii)=1
1755               l_homo(k,ii)=.true.
1756
1757 c             write (iout,*) "k",k
1758 c             write (iout,*) "i",i," j",j," constr_homology",
1759 c    &                       constr_homology
1760               ires_homo(ii)=i
1761               jres_homo(ii)=j
1762               odl(k,ii)=distal
1763               if (read2sigma) then
1764                 sigma_odl(k,ii)=0
1765                 do ik=i,j
1766                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1767                 enddo
1768                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1769                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1770      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1771               else
1772                 if (odl(k,ii).le.dist_cut) then
1773                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1774                 else
1775 #ifdef OLDSIGMA
1776                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1777      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1778 #else
1779                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1780      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1781 #endif
1782                 endif
1783               endif
1784               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1785             else
1786               ii=ii+1
1787 c              l_homo(k,ii)=.false.
1788             endif
1789             enddo
1790           enddo
1791         lim_odl=ii
1792         endif
1793 c
1794 c     Theta, dihedral and SC retraints
1795 c
1796         if (waga_angle.gt.0.0d0) then
1797           do i = nnt+3,nct
1798             if (idomain(k,i).eq.0) then
1799 c               sigma_dih(k,i)=0.0
1800                cycle
1801             endif
1802             dih(k,i)=phiref(i)
1803             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1804      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1805 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1806 c     &       " sigma_dihed",sigma_dih(k,i)
1807             if (sigma_dih(k,i).ne.0)
1808      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1809           enddo
1810           lim_dih=nct-nnt-2
1811         endif
1812
1813         if (waga_theta.gt.0.0d0) then
1814           do i = nnt+2,nct
1815              if (idomain(k,i).eq.0) then
1816 c              sigma_theta(k,i)=0.0
1817               cycle
1818              endif
1819              thetatpl(k,i)=thetaref(i)
1820              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1821      &                        rescore(k,i-2))/3.0
1822              if (sigma_theta(k,i).ne.0)
1823      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1824           enddo
1825         endif
1826
1827         if (waga_d.gt.0.0d0) then
1828           do i = nnt,nct
1829                if (itype(i).eq.10) cycle
1830                if (idomain(k,i).eq.0 ) then
1831 c                  sigma_d(k,i)=0.0
1832                   cycle
1833                endif
1834                xxtpl(k,i)=xxref(i)
1835                yytpl(k,i)=yyref(i)
1836                zztpl(k,i)=zzref(i)
1837                sigma_d(k,i)=rescore(k,i)
1838                if (sigma_d(k,i).ne.0)
1839      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1840                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1841           enddo
1842         endif
1843       enddo ! l
1844       enddo ! ll
1845 c
1846 c remove distance restraints not used in any model from the list
1847 c shift data in all arrays
1848 c
1849       if (waga_dist.ne.0.0d0) then
1850         ii=0
1851         liiflag=.true.
1852         do i=nnt,nct-2
1853          do j=i+2,nct
1854           ii=ii+1
1855           if (ii_in_use(ii).eq.0.and.liiflag) then
1856             liiflag=.false.
1857             iistart=ii
1858           endif
1859           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1860      &                   .not.liiflag.and.ii.eq.lim_odl) then
1861              if (ii.eq.lim_odl) then
1862               iishift=ii-iistart+1
1863              else
1864               iishift=ii-iistart
1865              endif
1866              liiflag=.true.
1867              do ki=iistart,lim_odl-iishift
1868               ires_homo(ki)=ires_homo(ki+iishift)
1869               jres_homo(ki)=jres_homo(ki+iishift)
1870               ii_in_use(ki)=ii_in_use(ki+iishift)
1871               do k=1,constr_homology
1872                odl(k,ki)=odl(k,ki+iishift)
1873                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1874                l_homo(k,ki)=l_homo(k,ki+iishift)
1875               enddo
1876              enddo
1877              ii=ii-iishift
1878              lim_odl=lim_odl-iishift
1879           endif
1880          enddo
1881         enddo
1882       endif
1883
1884       return
1885    10 stop "Error infragment file"
1886       end
1887 c----------------------------------------------------------------------