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