source/cluster/wham/src-M-homology/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       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
  31       call readi(controlcard,'NRES',nres,0)
  32       call readi(controlcard,'RESCALE',rescale_mode,2)
  33       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
  34       write (iout,*) "DISTCHAINMAX",distchainmax
  35 C Reading the dimensions of box in x,y,z coordinates
  36       call reada(controlcard,'BOXX',boxxsize,100.0d0)
  37       call reada(controlcard,'BOXY',boxysize,100.0d0)
  38       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
  39 c Cutoff range for interactions
  40       call reada(controlcard,"R_CUT",r_cut,15.0d0)
  41       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
  42       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
  43       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
  44       if (lipthick.gt.0.0d0) then
  45        bordliptop=(boxzsize+lipthick)/2.0
  46        bordlipbot=bordliptop-lipthick
  47 C      endif
  48       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
  49      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
  50       buflipbot=bordlipbot+lipbufthick
  51       bufliptop=bordliptop-lipbufthick
  52       if ((lipbufthick*2.0d0).gt.lipthick)
  53      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
  54       endif
  55       write(iout,*) "bordliptop=",bordliptop
  56       write(iout,*) "bordlipbot=",bordlipbot
  57       write(iout,*) "bufliptop=",bufliptop
  58       write(iout,*) "buflipbot=",buflipbot
  59 C Shielding mode
  60       call readi(controlcard,'SHIELD',shield_mode,0)
  61       write (iout,*) "SHIELD MODE",shield_mode
  62       if (shield_mode.gt.0) then
  63       pi=3.141592d0
  64 C VSolvSphere the volume of solving sphere
  65 C      print *,pi,"pi"
  66 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
  67 C there will be no distinction between proline peptide group and normal peptide
  68 C group in case of shielding parameters
  69       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
  70       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
  71       write (iout,*) VSolvSphere,VSolvSphere_div
  72 C long axis of side chain
  73       do i=1,ntyp
  74       long_r_sidechain(i)=vbldsc0(1,i)
  75       short_r_sidechain(i)=sigma0(i)
  76       enddo
  77       buff_shield=1.0d0
  78       endif
  79       call readi(controlcard,'PDBOUT',outpdb,0)
  80       call readi(controlcard,'MOL2OUT',outmol2,0)
  81       refstr=(index(controlcard,'REFSTR').gt.0)
  82       write (iout,*) "REFSTR",refstr
  83       pdbref=(index(controlcard,'PDBREF').gt.0)
  84       iscode=index(controlcard,'ONE_LETTER')
  85       tree=(index(controlcard,'MAKE_TREE').gt.0)
  86       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  87       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  88       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  89      & " CONSTR_DIST",constr_dist
  90       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
  91       write (iout,*) "with_theta_constr ",with_theta_constr
  92       call flush(iout)
  93       min_var=(index(controlcard,'MINVAR').gt.0)
  94       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
  95       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
  96       print_fittest=(index(controlcard,'PRINT_FITTEST').gt.0)
  97       call readi(controlcard,'NCUT',ncut,0)
  98       call readi(controlcard,'NCLUST',nclust,5)
  99       call readi(controlcard,'SYM',symetr,1)
 100       write (iout,*) 'sym', symetr
 101       call readi(controlcard,'NSTART',nstart,0)
 102       call readi(controlcard,'NEND',nend,0)
 103       call reada(controlcard,'ECUT',ecut,10.0d0)
 104       call reada(controlcard,'PROB',prob_limit,0.99d0)
 105       write (iout,*) "Probability limit",prob_limit
 106       lgrp=(index(controlcard,'LGRP').gt.0)
 107       caonly=(index(controlcard,'CA_ONLY').gt.0)
 108       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
 109       if (ncut.gt.0)
 110      & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
 111       call readi(controlcard,'IOPT',iopt,2)
 112       lside = index(controlcard,"SIDE").gt.0
 113       lefree = index(controlcard,"EFREE").gt.0
 114       call readi(controlcard,'NTEMP',nT,1)
 115       write (iout,*) "nT",nT
 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       if (restr_type.eq.12)
 881      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
 882       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 883       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 884       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 885       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 886       normalize = index(controlcard,"NORMALIZE").gt.0
 887       write (iout,*) "WBOLTZD",wboltzd
 888 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 889 c      write (iout,*) "IFRAG"
 890 c      do i=1,nfrag_
 891 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 892 c      enddo
 893 c      write (iout,*) "IPAIR"
 894 c      do i=1,npair_
 895 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 896 c      enddo
 897       if (nfrag_.gt.0) then
 898         nres0=nres
 899         read(inp,'(a)') pdbfile
 900         write (iout,*)
 901      & "Distance restraints will be constructed from structure ",pdbfile
 902         open(ipdbin,file=pdbfile,status='old',err=11)
 903         call readpdb(.true.)
 904         nres=nres0
 905         close(ipdbin)
 906       endif
 907       do i=1,nfrag_
 908         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 909         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 910      &    ifrag_(2,i)=nstart_sup+nsup-1
 911 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 912 c        call flush(iout)
 913         if (wfrag_(i).eq.0.0d0) cycle
 914         do j=ifrag_(1,i),ifrag_(2,i)-1
 915           do k=j+1,ifrag_(2,i)
 916 c            write (iout,*) "j",j," k",k
 917             ddjk=dist(j,k)
 918             if (restr_type.eq.1) then
 919               nhpb=nhpb+1
 920               irestr_type(nhpb)=1
 921               ihpb(nhpb)=j
 922               jhpb(nhpb)=k
 923               dhpb(nhpb)=ddjk
 924               forcon(nhpb)=wfrag_(i)
 925             else if (constr_dist.eq.2) then
 926               if (ddjk.le.dist_cut) then
 927                 nhpb=nhpb+1
 928                 irestr_type(nhpb)=1
 929                 ihpb(nhpb)=j
 930                 jhpb(nhpb)=k
 931                 dhpb(nhpb)=ddjk
 932                 forcon(nhpb)=wfrag_(i)
 933               endif
 934             else if (restr_type.eq.3) then
 935               nhpb=nhpb+1
 936               irestr_type(nhpb)=1
 937               ihpb(nhpb)=j
 938               jhpb(nhpb)=k
 939               dhpb(nhpb)=ddjk
 940               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 941             endif
 942             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
 943      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 944           enddo
 945         enddo
 946       enddo
 947       do i=1,npair_
 948         if (wpair_(i).eq.0.0d0) cycle
 949         ii = ipair_(1,i)
 950         jj = ipair_(2,i)
 951         if (ii.gt.jj) then
 952           itemp=ii
 953           ii=jj
 954           jj=itemp
 955         endif
 956         do j=ifrag_(1,ii),ifrag_(2,ii)
 957           do k=ifrag_(1,jj),ifrag_(2,jj)
 958             if (restr_type.eq.1) then
 959               nhpb=nhpb+1
 960               irestr_type(nhpb)=1
 961               ihpb(nhpb)=j
 962               jhpb(nhpb)=k
 963               dhpb(nhpb)=ddjk
 964               forcon(nhpb)=wfrag_(i)
 965             else if (constr_dist.eq.2) then
 966               if (ddjk.le.dist_cut) then
 967                 nhpb=nhpb+1
 968                 irestr_type(nhpb)=1
 969                 ihpb(nhpb)=j
 970                 jhpb(nhpb)=k
 971                 dhpb(nhpb)=ddjk
 972                 forcon(nhpb)=wfrag_(i)
 973               endif
 974             else if (restr_type.eq.3) then
 975               nhpb=nhpb+1
 976               irestr_type(nhpb)=1
 977               ihpb(nhpb)=j
 978               jhpb(nhpb)=k
 979               dhpb(nhpb)=ddjk
 980               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 981             endif
 982             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
 983      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 984           enddo
 985         enddo
 986       enddo
 987
 988 c      print *,ndist_
 989       write (iout,*) "Distance restraints as read from input"
 990       do i=1,ndist_
 991         if (restr_type.eq.12) then
 992           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
 993      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
 994      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
 995      &    fordepth_peak(nhpb_peak+1),npeak
 996 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
 997 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
 998 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
 999 c     &    fordepth_peak(nhpb_peak+1),npeak
1000           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1001      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1002           nhpb_peak=nhpb_peak+1
1003           irestr_type_peak(nhpb_peak)=12
1004           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1005           ipeak(2,npeak)=i
1006           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1007      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1008      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1009      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1010      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1011           if (ibecarb_peak(nhpb_peak).gt.0) then
1012             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1013             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1014           endif
1015         else if (restr_type.eq.11) then
1016           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1017      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1018 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1019           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1020           nhpb=nhpb+1
1021           irestr_type(nhpb)=11
1022           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1023      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1024      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1025           if (ibecarb(nhpb).gt.0) then
1026             ihpb(nhpb)=ihpb(nhpb)+nres
1027             jhpb(nhpb)=jhpb(nhpb)+nres
1028           endif
1029         else if (constr_dist.eq.10) then
1030 c Cross-lonk Markov-like potential
1031           call card_concat(controlcard)
1032           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1033           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1034           ibecarb(nhpb+1)=0
1035           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1036           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1037           if (index(controlcard,"ZL").gt.0) then
1038             link_type=1
1039           else if (index(controlcard,"ADH").gt.0) then
1040             link_type=2
1041           else if (index(controlcard,"PDH").gt.0) then
1042             link_type=3
1043           else if (index(controlcard,"DSS").gt.0) then
1044             link_type=4
1045           else
1046             link_type=0
1047           endif
1048           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1049      &       xlink(1,link_type))
1050           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1051      &       xlink(2,link_type))
1052           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1053      &       xlink(3,link_type))
1054           call reada(controlcard,"SIGMA",forcon(nhpb+1),
1055      &       xlink(4,link_type))
1056           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1057 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1058 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1059           if (forcon(nhpb+1).le.0.0d0 .or.
1060      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1061           nhpb=nhpb+1
1062           irestr_type(nhpb)=10
1063           if (ibecarb(nhpb).gt.0) then
1064             ihpb(nhpb)=ihpb(nhpb)+nres
1065             jhpb(nhpb)=jhpb(nhpb)+nres
1066           endif
1067           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1068      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1069      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1070      &     irestr_type(nhpb)
1071         else
1072 C        print *,"in else"
1073           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1074      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1075           if (forcon(nhpb+1).gt.0.0d0) then
1076           nhpb=nhpb+1
1077           if (dhpb1(nhpb).eq.0.0d0) then
1078             irestr_type(nhpb)=1
1079           else
1080             irestr_type(nhpb)=2
1081           endif
1082           if (ibecarb(nhpb).gt.0) then
1083             ihpb(nhpb)=ihpb(nhpb)+nres
1084             jhpb(nhpb)=jhpb(nhpb)+nres
1085           endif
1086           if (dhpb(nhpb).eq.0.0d0)
1087      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1088         endif
1089         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1090      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1091         endif
1092 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1093 C        if (forcon(nhpb+1).gt.0.0d0) then
1094 C          nhpb=nhpb+1
1095 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1096       enddo
1097
1098       ENDDO ! next
1099
1100       fordepthmax=0.0d0
1101       if (normalize) then
1102         do i=nss+1,nhpb
1103           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1104      &      fordepthmax=fordepth(i)
1105         enddo
1106         do i=nss+1,nhpb
1107           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1108         enddo
1109       endif
1110       if (nhpb.gt.nss)  then
1111         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1112      &  "The following",nhpb-nss,
1113      &  " distance restraints have been imposed:",
1114      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
1115      &  "  score"," type"
1116         do i=nss+1,nhpb
1117           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1118      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1119      &  irestr_type(i)
1120         enddo
1121       endif
1122       call hpb_partition
1123       call flush(iout)
1124       return
1125    11 write (iout,*)"read_dist_restr: error reading reference structure"
1126       stop
1127       end
1128 c-------------------------------------------------------------------------------
1129       subroutine read_saxs_constr
1130       implicit real*8 (a-h,o-z)
1131       include 'DIMENSIONS'
1132 #ifdef MPI
1133       include 'mpif.h'
1134 #endif
1135       include 'COMMON.CONTROL'
1136       include 'COMMON.CHAIN'
1137       include 'COMMON.IOUNITS'
1138       include 'COMMON.SBRIDGE'
1139       double precision cm(3)
1140 c      read(inp,*) nsaxs
1141       write (iout,*) "Calling read_saxs nsaxs",nsaxs
1142       call flush(iout)
1143       if (saxs_mode.eq.0) then
1144 c SAXS distance distribution
1145       do i=1,nsaxs
1146         read(inp,*) distsaxs(i),Psaxs(i)
1147       enddo
1148       Cnorm = 0.0d0
1149       do i=1,nsaxs
1150         Cnorm = Cnorm + Psaxs(i)
1151       enddo
1152       write (iout,*) "Cnorm",Cnorm
1153       do i=1,nsaxs
1154         Psaxs(i)=Psaxs(i)/Cnorm
1155       enddo
1156       write (iout,*) "Normalized distance distribution from SAXS"
1157       do i=1,nsaxs
1158         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1159       enddo
1160       Wsaxs0=0.0d0
1161       do i=1,nsaxs
1162         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1163       enddo
1164       write (iout,*) "Wsaxs0",Wsaxs0
1165       else
1166 c SAXS "spheres".
1167       do i=1,nsaxs
1168         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1169       enddo
1170       do j=1,3
1171         cm(j)=0.0d0
1172       enddo
1173       do i=1,nsaxs
1174         do j=1,3
1175           cm(j)=cm(j)+Csaxs(j,i)
1176         enddo
1177       enddo
1178       do j=1,3
1179         cm(j)=cm(j)/nsaxs
1180       enddo
1181       do i=1,nsaxs
1182         do j=1,3
1183           Csaxs(j,i)=Csaxs(j,i)-cm(j)
1184         enddo
1185       enddo
1186       write (iout,*) "SAXS sphere coordinates"
1187       do i=1,nsaxs
1188         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
1189       enddo
1190       endif
1191       return
1192       end
1193 c====-------------------------------------------------------------------
1194       subroutine read_constr_homology(lprn)
1195
1196       include 'DIMENSIONS'
1197 #ifdef MPI
1198       include 'mpif.h'
1199 #endif
1200       include 'COMMON.SETUP'
1201       include 'COMMON.CONTROL'
1202       include 'COMMON.CHAIN'
1203       include 'COMMON.IOUNITS'
1204       include 'COMMON.GEO'
1205       include 'COMMON.INTERACT'
1206       include 'COMMON.NAMES'
1207       include 'COMMON.HOMRESTR'
1208 c
1209 c For new homol impl
1210 c
1211       include 'COMMON.VAR'
1212 c     include 'include_unres/COMMON.VAR'
1213 c
1214
1215 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1216 c    &                 dist_cut
1217 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1218 c    &    sigma_odl_temp(maxres,maxres,max_template)
1219       character*2 kic2
1220       character*24 model_ki_dist, model_ki_angle
1221       character*500 controlcard
1222       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1223       integer idomain(max_template,maxres)
1224       integer ilen
1225       external ilen
1226       logical lprn
1227       logical unres_pdb,liiflag
1228 c
1229 c     FP - Nov. 2014 Temporary specifications for new vars
1230 c
1231       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
1232      &                 rescore3_tmp
1233       double precision, dimension (max_template,maxres) :: rescore
1234       double precision, dimension (max_template,maxres) :: rescore2
1235       double precision, dimension (max_template,maxres) :: rescore3
1236       character*24 tpl_k_rescore
1237 c -----------------------------------------------------------------
1238 c Reading multiple PDB ref structures and calculation of retraints
1239 c not using pre-computed ones stored in files model_ki_{dist,angle}
1240 c FP (Nov., 2014)
1241 c -----------------------------------------------------------------
1242 c
1243 c
1244 c Alternative: reading from input
1245 #ifdef DEBUG
1246       write (iout,*) "BEGIN READ HOMOLOGY INFO"
1247 #ifdef AIX
1248       call flush_(iout)
1249 #else
1250       call flush(iout)
1251 #endif
1252 #endif
1253       call card_concat(controlcard)
1254       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1255       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1256       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1257       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1258       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1259       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1260       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1261       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1262       if (homol_nset.gt.1)then
1263          call readi(controlcard,"ISET",iset,1)
1264          call card_concat(controlcard)
1265          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1266       else
1267         iset=1
1268         waga_homology(1)=1.0
1269       endif
1270 c
1271 #ifdef DEBUG
1272       write(iout,*) "read_constr_homology iset",iset
1273       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1274 #ifdef AIX
1275       call flush_(iout)
1276 #else
1277       call flush(iout)
1278 #endif
1279 #endif
1280 cd      write (iout,*) "nnt",nnt," nct",nct
1281 cd      call flush(iout)
1282
1283
1284       lim_odl=0
1285       lim_dih=0
1286 c
1287 c  New
1288 c
1289 c
1290 c  Reading HM global scores (prob not required)
1291 c
1292       do i = nnt,nct
1293         do k=1,constr_homology
1294          idomain(k,i)=0
1295         enddo
1296       enddo
1297 c     open (4,file="HMscore")
1298 c     do k=1,constr_homology
1299 c       read (4,*,end=521) hmscore_tmp
1300 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1301 c       write(*,*) "Model", k, ":", hmscore(k)
1302 c     enddo
1303 c521  continue
1304
1305       ii=0
1306       do i = nnt,nct-2
1307         do j=i+2,nct
1308         ii=ii+1
1309         ii_in_use(ii)=0
1310         enddo
1311       enddo
1312 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1313
1314       if (read_homol_frag) then
1315         call read_klapaucjusz
1316       else
1317       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1318       do k=1,constr_homology
1319
1320         read(inp,'(a)') pdbfile
1321 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1322 c  Next stament causes error upon compilation (?)
1323 c       if(me.eq.king.or. .not. out1file)
1324 c         write (iout,'(2a)') 'PDB data will be read from file ',
1325 c    &   pdbfile(:ilen(pdbfile))
1326          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1327      &  pdbfile(:ilen(pdbfile))
1328         open(ipdbin,file=pdbfile,status='old',err=33)
1329         goto 34
1330   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1331      &  pdbfile(:ilen(pdbfile))
1332         stop
1333   34    continue
1334 c        print *,'Begin reading pdb data'
1335 c
1336 c Files containing res sim or local scores (former containing sigmas)
1337 c
1338
1339         write(kic2,'(bz,i2.2)') k
1340
1341         tpl_k_rescore="template"//kic2//".sco"
1342
1343         unres_pdb=.false.
1344         call readpdb_template(k)
1345         do i=1,2*nres
1346           do j=1,3
1347             crefjlee(j,i)=c(j,i)
1348           enddo
1349         enddo
1350 #ifdef DEBUG
1351         do i=1,nres
1352           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1353      &      (crefjlee(j,i+nres),j=1,3)
1354         enddo
1355         write (iout,*) "READ HOMOLOGY INFO"
1356         write (iout,*) "read_constr_homology x: after reading pdb file"
1357         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1358         write (iout,*) "waga_dist",waga_dist
1359         write (iout,*) "waga_angle",waga_angle
1360         write (iout,*) "waga_theta",waga_theta
1361         write (iout,*) "waga_d",waga_d
1362         write (iout,*) "dist_cut",dist_cut
1363 #endif
1364 #ifdef AIX
1365       call flush_(iout)
1366 #else
1367       call flush(iout)
1368 #endif
1369
1370 c
1371 c     Distance restraints
1372 c
1373 c          ... --> odl(k,ii)
1374 C Copy the coordinates from reference coordinates (?)
1375         do i=1,2*nres
1376           do j=1,3
1377 c            c(j,i)=cref(j,i)
1378 c           write (iout,*) "c(",j,i,") =",c(j,i)
1379           enddo
1380         enddo
1381 c
1382 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1383 c
1384 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1385           open (ientin,file=tpl_k_rescore,status='old')
1386           if (nnt.gt.1) rescore(k,1)=0.0d0
1387           do irec=nnt,nct ! loop for reading res sim
1388             if (read2sigma) then
1389              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1390      &                                rescore3_tmp,idomain_tmp
1391              i_tmp=i_tmp+nnt-1
1392              idomain(k,i_tmp)=idomain_tmp
1393              rescore(k,i_tmp)=rescore_tmp
1394              rescore2(k,i_tmp)=rescore2_tmp
1395              rescore3(k,i_tmp)=rescore3_tmp
1396              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1397      &                      i_tmp,rescore2_tmp,rescore_tmp,
1398      &                                rescore3_tmp,idomain_tmp
1399             else
1400              idomain(k,irec)=1
1401              read (ientin,*,end=1401) rescore_tmp
1402
1403 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1404              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1405 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1406             endif
1407           enddo
1408  1401   continue
1409         close (ientin)
1410         if (waga_dist.ne.0.0d0) then
1411           ii=0
1412           do i = nnt,nct-2
1413             do j=i+2,nct
1414
1415               x12=c(1,i)-c(1,j)
1416               y12=c(2,i)-c(2,j)
1417               z12=c(3,i)-c(3,j)
1418               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1419 c              write (iout,*) k,i,j,distal,dist2_cut
1420
1421             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1422      &            .and. distal.le.dist2_cut ) then
1423
1424               ii=ii+1
1425               ii_in_use(ii)=1
1426               l_homo(k,ii)=.true.
1427
1428 c             write (iout,*) "k",k
1429 c             write (iout,*) "i",i," j",j," constr_homology",
1430 c    &                       constr_homology
1431               ires_homo(ii)=i
1432               jres_homo(ii)=j
1433               odl(k,ii)=distal
1434               if (read2sigma) then
1435                 sigma_odl(k,ii)=0
1436                 do ik=i,j
1437                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1438                 enddo
1439                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1440                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1441      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1442               else
1443                 if (odl(k,ii).le.dist_cut) then
1444                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1445                 else
1446 #ifdef OLDSIGMA
1447                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1448      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1449 #else
1450                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1451      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1452 #endif
1453                 endif
1454               endif
1455               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1456             else
1457               ii=ii+1
1458               l_homo(k,ii)=.false.
1459             endif
1460             enddo
1461           enddo
1462         lim_odl=ii
1463         endif
1464 c
1465 c     Theta, dihedral and SC retraints
1466 c
1467         if (waga_angle.gt.0.0d0) then
1468 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1469 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1470 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1471 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1472 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1473 c    &                            sigma_dih(k,i+nnt-1)
1474 c         enddo
1475 c1402   continue
1476 c         close (ientin)
1477           do i = nnt+3,nct
1478             if (idomain(k,i).eq.0) then
1479                sigma_dih(k,i)=0.0
1480                cycle
1481             endif
1482             dih(k,i)=phiref(i) ! right?
1483 c           read (ientin,*) sigma_dih(k,i) ! original variant
1484 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1485 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1486 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1487 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1488 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1489
1490             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1491      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1492 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1493 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1494 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1495 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1496 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1497 c   Instead of res sim other local measure of b/b str reliability possible
1498             if (sigma_dih(k,i).ne.0)
1499      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1500 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1501           enddo
1502           lim_dih=nct-nnt-2
1503         endif
1504
1505         if (waga_theta.gt.0.0d0) then
1506 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1507 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1508 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1509 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1510 c    &                              sigma_theta(k,i+nnt-1)
1511 c         enddo
1512 c1403   continue
1513 c         close (ientin)
1514
1515           do i = nnt+2,nct ! right? without parallel.
1516 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1517 c         do i=ithet_start,ithet_end ! with FG parallel.
1518              if (idomain(k,i).eq.0) then
1519               sigma_theta(k,i)=0.0
1520               cycle
1521              endif
1522              thetatpl(k,i)=thetaref(i)
1523 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1524 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1525 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1526 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1527 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1528              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1529      &                        rescore(k,i-2))/3.0
1530 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1531              if (sigma_theta(k,i).ne.0)
1532      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1533
1534 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1535 c                             rescore(k,i-2) !  right expression ?
1536 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1537           enddo
1538         endif
1539
1540         if (waga_d.gt.0.0d0) then
1541 c       open (ientin,file=tpl_k_sigma_d,status='old')
1542 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1543 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1544 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1545 c    &                          sigma_d(k,i+nnt-1)
1546 c         enddo
1547 c1404   continue
1548
1549           do i = nnt,nct ! right? without parallel.
1550 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1551 c         do i=loc_start,loc_end ! with FG parallel.
1552                if (itype(i).eq.10) cycle
1553                if (idomain(k,i).eq.0 ) then
1554                   sigma_d(k,i)=0.0
1555                   cycle
1556                endif
1557                xxtpl(k,i)=xxref(i)
1558                yytpl(k,i)=yyref(i)
1559                zztpl(k,i)=zzref(i)
1560 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1561 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1562 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1563 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1564                sigma_d(k,i)=rescore3(k,i) !  right expression ?
1565                if (sigma_d(k,i).ne.0)
1566      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1567
1568 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1569 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1570 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1571           enddo
1572         endif
1573       enddo
1574 c
1575 c remove distance restraints not used in any model from the list
1576 c shift data in all arrays
1577 c
1578       if (waga_dist.ne.0.0d0) then
1579         ii=0
1580         liiflag=.true.
1581         do i=nnt,nct-2
1582          do j=i+2,nct
1583           ii=ii+1
1584           if (ii_in_use(ii).eq.0.and.liiflag) then
1585             liiflag=.false.
1586             iistart=ii
1587           endif
1588           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1589      &                   .not.liiflag.and.ii.eq.lim_odl) then
1590              if (ii.eq.lim_odl) then
1591               iishift=ii-iistart+1
1592              else
1593               iishift=ii-iistart
1594              endif
1595              liiflag=.true.
1596              do ki=iistart,lim_odl-iishift
1597               ires_homo(ki)=ires_homo(ki+iishift)
1598               jres_homo(ki)=jres_homo(ki+iishift)
1599               ii_in_use(ki)=ii_in_use(ki+iishift)
1600               do k=1,constr_homology
1601                odl(k,ki)=odl(k,ki+iishift)
1602                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1603                l_homo(k,ki)=l_homo(k,ki+iishift)
1604               enddo
1605              enddo
1606              ii=ii-iishift
1607              lim_odl=lim_odl-iishift
1608           endif
1609          enddo
1610         enddo
1611       endif
1612
1613       endif ! .not. klapaucjusz
1614
1615       if (constr_homology.gt.0) call homology_partition
1616       if (constr_homology.gt.0) call init_int_table
1617 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1618 cd     & "lim_xx=",lim_xx
1619 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1620 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1621 c
1622 c Print restraints
1623 c
1624       if (.not.lprn) return
1625 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1626       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1627        write (iout,*) "Distance restraints from templates"
1628        do ii=1,lim_odl
1629        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1630      &  ii,ires_homo(ii),jres_homo(ii),
1631      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1632      &  ki=1,constr_homology)
1633        enddo
1634        write (iout,*) "Dihedral angle restraints from templates"
1635        do i=nnt+3,nct
1636         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1637      &      (rad2deg*dih(ki,i),
1638      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1639        enddo
1640        write (iout,*) "Virtual-bond angle restraints from templates"
1641        do i=nnt+2,nct
1642         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1643      &      (rad2deg*thetatpl(ki,i),
1644      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1645        enddo
1646        write (iout,*) "SC restraints from templates"
1647        do i=nnt,nct
1648         write(iout,'(i5,100(4f8.2,4x))') i,
1649      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1650      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1651        enddo
1652       endif
1653 c -----------------------------------------------------------------
1654       return
1655       end
1656 c----------------------------------------------------------------------
1657       subroutine read_klapaucjusz
1658
1659       include 'DIMENSIONS'
1660 #ifdef MPI
1661       include 'mpif.h'
1662 #endif
1663       include 'COMMON.SETUP'
1664       include 'COMMON.CONTROL'
1665       include 'COMMON.CHAIN'
1666       include 'COMMON.IOUNITS'
1667       include 'COMMON.GEO'
1668       include 'COMMON.INTERACT'
1669       include 'COMMON.NAMES'
1670       include 'COMMON.HOMRESTR'
1671       character*256 fragfile
1672       integer ninclust(maxclust),inclust(max_template,maxclust),
1673      &  nresclust(maxclust),iresclust(maxres,maxclust)
1674
1675       character*2 kic2
1676       character*24 model_ki_dist, model_ki_angle
1677       character*500 controlcard
1678       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1679       integer idomain(max_template,maxres)
1680       logical lprn /.true./
1681       integer ilen
1682       external ilen
1683       logical unres_pdb,liiflag
1684 c
1685 c
1686       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1687       double precision, dimension (max_template,maxres) :: rescore
1688       double precision, dimension (max_template,maxres) :: rescore2
1689       character*24 tpl_k_rescore
1690
1691 c
1692 c For new homol impl
1693 c
1694       include 'COMMON.VAR'
1695 c
1696       double precision chomo(3,maxres2+2,max_template)
1697       call getenv("FRAGFILE",fragfile)
1698       open(ientin,file=fragfile,status="old",err=10)
1699       read(ientin,*) constr_homology,nclust
1700       l_homo = .false.
1701       sigma_theta=0.0
1702       sigma_d=0.0
1703       sigma_dih=0.0
1704 c Read pdb files
1705       do k=1,constr_homology
1706         read(ientin,'(a)') pdbfile
1707         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1708      &  pdbfile(:ilen(pdbfile))
1709         open(ipdbin,file=pdbfile,status='old',err=33)
1710         goto 34
1711   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1712      &  pdbfile(:ilen(pdbfile))
1713         stop
1714   34    continue
1715         unres_pdb=.false.
1716         call readpdb_template(k)
1717         do i=1,2*nres
1718           do j=1,3
1719             chomo(j,i,k)=c(j,i)
1720           enddo
1721         enddo
1722         do i=1,nres
1723           rescore(k,i)=0.2d0
1724           rescore2(k,i)=1.0d0
1725         enddo
1726       enddo
1727 c Read clusters
1728       do i=1,nclust
1729         read(ientin,*) ninclust(i),nresclust(i)
1730         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1731         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1732       enddo
1733 c
1734 c Loop over clusters
1735 c
1736       do l=1,nclust
1737         do ll = 1,ninclust(l)
1738
1739         k = inclust(ll,l)
1740         do i=1,nres
1741           idomain(k,i)=0
1742         enddo
1743         do i=1,nresclust(l)
1744           if (nnt.gt.1)  then
1745             idomain(k,iresclust(i,l)+1) = 1
1746           else
1747             idomain(k,iresclust(i,l)) = 1
1748           endif
1749         enddo
1750 c
1751 c     Distance restraints
1752 c
1753 c          ... --> odl(k,ii)
1754 C Copy the coordinates from reference coordinates (?)
1755         do i=1,2*nres
1756           do j=1,3
1757             c(j,i)=chomo(j,i,k)
1758 c           write (iout,*) "c(",j,i,") =",c(j,i)
1759           enddo
1760         enddo
1761         call int_from_cart1(.false.)
1762         call int_from_cart(.true.,.false.)
1763         call sc_loc_geom(.false.)
1764         do i=1,nres
1765           thetaref(i)=theta(i)
1766           phiref(i)=phi(i)
1767         enddo
1768         if (waga_dist.ne.0.0d0) then
1769           ii=0
1770           do i = nnt,nct-2
1771             do j=i+2,nct
1772
1773               x12=c(1,i)-c(1,j)
1774               y12=c(2,i)-c(2,j)
1775               z12=c(3,i)-c(3,j)
1776               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1777 c              write (iout,*) k,i,j,distal,dist2_cut
1778
1779             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1780      &            .and. distal.le.dist2_cut ) then
1781
1782               ii=ii+1
1783               ii_in_use(ii)=1
1784               l_homo(k,ii)=.true.
1785
1786 c             write (iout,*) "k",k
1787 c             write (iout,*) "i",i," j",j," constr_homology",
1788 c    &                       constr_homology
1789               ires_homo(ii)=i
1790               jres_homo(ii)=j
1791               odl(k,ii)=distal
1792               if (read2sigma) then
1793                 sigma_odl(k,ii)=0
1794                 do ik=i,j
1795                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1796                 enddo
1797                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1798                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1799      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1800               else
1801                 if (odl(k,ii).le.dist_cut) then
1802                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1803                 else
1804 #ifdef OLDSIGMA
1805                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1806      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1807 #else
1808                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1809      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1810 #endif
1811                 endif
1812               endif
1813               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1814             else
1815               ii=ii+1
1816 c              l_homo(k,ii)=.false.
1817             endif
1818             enddo
1819           enddo
1820         lim_odl=ii
1821         endif
1822 c
1823 c     Theta, dihedral and SC retraints
1824 c
1825         if (waga_angle.gt.0.0d0) then
1826           do i = nnt+3,nct
1827             if (idomain(k,i).eq.0) then
1828 c               sigma_dih(k,i)=0.0
1829                cycle
1830             endif
1831             dih(k,i)=phiref(i)
1832             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1833      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1834 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1835 c     &       " sigma_dihed",sigma_dih(k,i)
1836             if (sigma_dih(k,i).ne.0)
1837      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1838           enddo
1839           lim_dih=nct-nnt-2
1840         endif
1841
1842         if (waga_theta.gt.0.0d0) then
1843           do i = nnt+2,nct
1844              if (idomain(k,i).eq.0) then
1845 c              sigma_theta(k,i)=0.0
1846               cycle
1847              endif
1848              thetatpl(k,i)=thetaref(i)
1849              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1850      &                        rescore(k,i-2))/3.0
1851              if (sigma_theta(k,i).ne.0)
1852      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1853           enddo
1854         endif
1855
1856         if (waga_d.gt.0.0d0) then
1857           do i = nnt,nct
1858                if (itype(i).eq.10) cycle
1859                if (idomain(k,i).eq.0 ) then
1860 c                  sigma_d(k,i)=0.0
1861                   cycle
1862                endif
1863                xxtpl(k,i)=xxref(i)
1864                yytpl(k,i)=yyref(i)
1865                zztpl(k,i)=zzref(i)
1866                sigma_d(k,i)=rescore(k,i)
1867                if (sigma_d(k,i).ne.0)
1868      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1869                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1870           enddo
1871         endif
1872       enddo ! l
1873       enddo ! ll
1874 c
1875 c remove distance restraints not used in any model from the list
1876 c shift data in all arrays
1877 c
1878       if (waga_dist.ne.0.0d0) then
1879         ii=0
1880         liiflag=.true.
1881         do i=nnt,nct-2
1882          do j=i+2,nct
1883           ii=ii+1
1884           if (ii_in_use(ii).eq.0.and.liiflag) then
1885             liiflag=.false.
1886             iistart=ii
1887           endif
1888           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1889      &                   .not.liiflag.and.ii.eq.lim_odl) then
1890              if (ii.eq.lim_odl) then
1891               iishift=ii-iistart+1
1892              else
1893               iishift=ii-iistart
1894              endif
1895              liiflag=.true.
1896              do ki=iistart,lim_odl-iishift
1897               ires_homo(ki)=ires_homo(ki+iishift)
1898               jres_homo(ki)=jres_homo(ki+iishift)
1899               ii_in_use(ki)=ii_in_use(ki+iishift)
1900               do k=1,constr_homology
1901                odl(k,ki)=odl(k,ki+iishift)
1902                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1903                l_homo(k,ki)=l_homo(k,ki+iishift)
1904               enddo
1905              enddo
1906              ii=ii-iishift
1907              lim_odl=lim_odl-iishift
1908           endif
1909          enddo
1910         enddo
1911       endif
1912
1913       return
1914    10 stop "Error infragment file"
1915       end
1916 c----------------------------------------------------------------------