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.FREE'
  17       include 'COMMON.INTERACT'
  18       include "COMMON.SPLITELE"
  19       include 'COMMON.SHIELD'
  20       character*320 controlcard,ucase
  21 #ifdef MPL
  22       include 'COMMON.INFO'
  23 #endif
  24       integer i,i1,i2,it1,it2
  25       double precision pi
  26       read (INP,'(a80)') titel
  27       call card_concat(controlcard)
  28
  29       call readi(controlcard,'NRES',nres,0)
  30       call readi(controlcard,'RESCALE',rescale_mode,2)
  31       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
  32       write (iout,*) "DISTCHAINMAX",distchainmax
  33 C Reading the dimensions of box in x,y,z coordinates
  34       call reada(controlcard,'BOXX',boxxsize,100.0d0)
  35       call reada(controlcard,'BOXY',boxysize,100.0d0)
  36       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
  37 c Cutoff range for interactions
  38       call reada(controlcard,"R_CUT",r_cut,15.0d0)
  39       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
  40       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
  41       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
  42       if (lipthick.gt.0.0d0) then
  43        bordliptop=(boxzsize+lipthick)/2.0
  44        bordlipbot=bordliptop-lipthick
  45 C      endif
  46       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
  47      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
  48       buflipbot=bordlipbot+lipbufthick
  49       bufliptop=bordliptop-lipbufthick
  50       if ((lipbufthick*2.0d0).gt.lipthick)
  51      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
  52       endif
  53       write(iout,*) "bordliptop=",bordliptop
  54       write(iout,*) "bordlipbot=",bordlipbot
  55       write(iout,*) "bufliptop=",bufliptop
  56       write(iout,*) "buflipbot=",buflipbot
  57 C Shielding mode
  58       call readi(controlcard,'SHIELD',shield_mode,0)
  59       write (iout,*) "SHIELD MODE",shield_mode
  60       if (shield_mode.gt.0) then
  61       pi=3.141592d0
  62 C VSolvSphere the volume of solving sphere
  63 C      print *,pi,"pi"
  64 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
  65 C there will be no distinction between proline peptide group and normal peptide
  66 C group in case of shielding parameters
  67       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
  68       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
  69       write (iout,*) VSolvSphere,VSolvSphere_div
  70 C long axis of side chain
  71       do i=1,ntyp
  72       long_r_sidechain(i)=vbldsc0(1,i)
  73       short_r_sidechain(i)=sigma0(i)
  74       enddo
  75       buff_shield=1.0d0
  76       endif
  77       call readi(controlcard,'PDBOUT',outpdb,0)
  78       call readi(controlcard,'MOL2OUT',outmol2,0)
  79       refstr=(index(controlcard,'REFSTR').gt.0)
  80       write (iout,*) "REFSTR",refstr
  81       pdbref=(index(controlcard,'PDBREF').gt.0)
  82       iscode=index(controlcard,'ONE_LETTER')
  83       tree=(index(controlcard,'MAKE_TREE').gt.0)
  84       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  85       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  86       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  87      & " CONSTR_DIST",constr_dist
  88       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
  89       write (iout,*) "with_theta_constr ",with_theta_constr
  90       call flush(iout)
  91       min_var=(index(controlcard,'MINVAR').gt.0)
  92       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
  93       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
  94       call readi(controlcard,'NCUT',ncut,1)
  95       call readi(controlcard,'SYM',symetr,1)
  96       write (iout,*) 'sym', symetr
  97       call readi(controlcard,'NSTART',nstart,0)
  98       call readi(controlcard,'NEND',nend,0)
  99       call reada(controlcard,'ECUT',ecut,10.0d0)
 100       call reada(controlcard,'PROB',prob_limit,0.99d0)
 101       write (iout,*) "Probability limit",prob_limit
 102       lgrp=(index(controlcard,'LGRP').gt.0)
 103       caonly=(index(controlcard,'CA_ONLY').gt.0)
 104       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
 105       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
 106       call readi(controlcard,'IOPT',iopt,2)
 107       lside = index(controlcard,"SIDE").gt.0
 108       efree = index(controlcard,"EFREE").gt.0
 109       call readi(controlcard,'NTEMP',nT,1)
 110       write (iout,*) "nT",nT
 111       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
 112       write (iout,*) "nT",nT
 113       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
 114       do i=1,nT
 115         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
 116       enddo
 117       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
 118       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
 119       lprint_int=index(controlcard,"PRINT_INT") .gt.0
 120       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
 121       write (iout,*) "with_homology_constr ",with_dihed_constr,
 122      & " CONSTR_HOMOLOGY",constr_homology
 123       print_homology_restraints=
 124      & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
 125       print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
 126       print_homology_models=
 127      & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
 128       if (min_var) iopt=1
 129       return
 130       end
 131 c--------------------------------------------------------------------------
 132       subroutine molread
 133 C
 134 C Read molecular data.
 135 C
 136       implicit none
 137       include 'DIMENSIONS'
 138       include 'COMMON.IOUNITS'
 139       include 'COMMON.GEO'
 140       include 'COMMON.VAR'
 141       include 'COMMON.INTERACT'
 142       include 'COMMON.LOCAL'
 143       include 'COMMON.NAMES'
 144       include 'COMMON.CHAIN'
 145       include 'COMMON.FFIELD'
 146       include 'COMMON.SBRIDGE'
 147       include 'COMMON.HEADER'
 148       include 'COMMON.CONTROL'
 149       include 'COMMON.CONTACTS'
 150       include 'COMMON.TIME1'
 151       include 'COMMON.TORCNSTR'
 152       include 'COMMON.SHIELD'
 153 #ifdef MPL
 154       include 'COMMON.INFO'
 155 #endif
 156       character*4 sequence(maxres)
 157       character*800 weightcard
 158       integer rescode
 159       double precision x(maxvar)
 160       integer itype_pdb(maxres)
 161       logical seq_comp
 162       integer i,j,kkk,i1,i2,it1,it2
 163 C
 164 C Body
 165 C
 166 C Read weights of the subsequent energy terms.
 167       call card_concat(weightcard)
 168       write(iout,*) weightcard
 169 C      call reada(weightcard,'WSC',wsc,1.0d0)
 170       write(iout,*) wsc
 171       call reada(weightcard,'WLONG',wsc,wsc)
 172       call reada(weightcard,'WSCP',wscp,1.0d0)
 173       call reada(weightcard,'WELEC',welec,1.0D0)
 174       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 175       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 176       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 177       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 178       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 179       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 180       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 181       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 182       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 183       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 184       call reada(weightcard,'WBOND',wbond,1.0D0)
 185       call reada(weightcard,'WTOR',wtor,1.0D0)
 186       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 187       call reada(weightcard,'WANG',wang,1.0D0)
 188       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 189       call reada(weightcard,'SCAL14',scal14,0.4D0)
 190       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 191       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 192       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 193       if (index(weightcard,'SOFT').gt.0) ipot=6
 194       call reada(weightcard,"D0CM",d0cm,3.78d0)
 195       call reada(weightcard,"AKCM",akcm,15.1d0)
 196       call reada(weightcard,"AKTH",akth,11.0d0)
 197       call reada(weightcard,"AKCT",akct,12.0d0)
 198       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 199       call reada(weightcard,"V2SS",v2ss,7.61d0)
 200       call reada(weightcard,"V3SS",v3ss,13.7d0)
 201       call reada(weightcard,"EBR",ebr,-5.50D0)
 202       call reada(weightcard,'WSHIELD',wshield,1.0d0)
 203       write(iout,*) 'WSHIELD',wshield
 204        call reada(weightcard,'WLT',wliptran,0.0D0)
 205       call reada(weightcard,"ATRISS",atriss,0.301D0)
 206       call reada(weightcard,"BTRISS",btriss,0.021D0)
 207       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 208       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 209       write (iout,*) "ATRISS=", atriss
 210       write (iout,*) "BTRISS=", btriss
 211       write (iout,*) "CTRISS=", ctriss
 212       write (iout,*) "DTRISS=", dtriss
 213       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 214       do i=1,maxres
 215         dyn_ss_mask(i)=.false.
 216       enddo
 217       do i=1,maxres-1
 218         do j=i+1,maxres
 219           dyn_ssbond_ij(i,j)=1.0d300
 220         enddo
 221       enddo
 222       call reada(weightcard,"HT",Ht,0.0D0)
 223       if (dyn_ss) then
 224         ss_depth=ebr/wsc-0.25*eps(1,1)
 225         Ht=Ht/wsc-0.25*eps(1,1)
 226         akcm=akcm*wstrain/wsc
 227         akth=akth*wstrain/wsc
 228         akct=akct*wstrain/wsc
 229         v1ss=v1ss*wstrain/wsc
 230         v2ss=v2ss*wstrain/wsc
 231         v3ss=v3ss*wstrain/wsc
 232       else
 233         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 234       endif
 235       write (iout,'(/a)') "Disulfide bridge parameters:"
 236       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 237       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 238       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 239       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 240      & ' v3ss:',v3ss
 241
 242 C 12/1/95 Added weight for the multi-body term WCORR
 243       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 244       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 245       weights(1)=wsc
 246       weights(2)=wscp
 247       weights(3)=welec
 248       weights(4)=wcorr
 249       weights(5)=wcorr5
 250       weights(6)=wcorr6
 251       weights(7)=wel_loc
 252       weights(8)=wturn3
 253       weights(9)=wturn4
 254       weights(10)=wturn6
 255       weights(11)=wang
 256       weights(12)=wscloc
 257       weights(13)=wtor
 258       weights(14)=wtor_d
 259       weights(15)=wstrain
 260       weights(16)=wvdwpp
 261       weights(17)=wbond
 262       weights(18)=scal14
 263       weights(19)=wsccor
 264       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 265      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 266      &  wturn4,wturn6,wsccor
 267    10 format (/'Energy-term weights (unscaled):'//
 268      & 'WSCC=   ',f10.6,' (SC-SC)'/
 269      & 'WSCP=   ',f10.6,' (SC-p)'/
 270      & 'WELEC=  ',f10.6,' (p-p electr)'/
 271      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 272      & 'WBOND=  ',f10.6,' (stretching)'/
 273      & 'WANG=   ',f10.6,' (bending)'/
 274      & 'WSCLOC= ',f10.6,' (SC local)'/
 275      & 'WTOR=   ',f10.6,' (torsional)'/
 276      & 'WTORD=  ',f10.6,' (double torsional)'/
 277      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 278      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 279      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 280      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 281      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 282      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 283      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 284      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 285      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 286
 287       if (wcorr4.gt.0.0d0) then
 288         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 289      &   'between contact pairs of peptide groups'
 290         write (iout,'(2(a,f5.3/))')
 291      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 292      &  'Range of quenching the correlation terms:',2*delt_corr
 293       else if (wcorr.gt.0.0d0) then
 294         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 295      &   'between contact pairs of peptide groups'
 296       endif
 297       write (iout,'(a,f8.3)')
 298      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 299       write (iout,'(a,f8.3)')
 300      &  'General scaling factor of SC-p interactions:',scalscp
 301       r0_corr=cutoff_corr-delt_corr
 302       do i=1,20
 303         aad(i,1)=scalscp*aad(i,1)
 304         aad(i,2)=scalscp*aad(i,2)
 305         bad(i,1)=scalscp*bad(i,1)
 306         bad(i,2)=scalscp*bad(i,2)
 307       enddo
 308
 309       call flush(iout)
 310       print *,'indpdb=',indpdb,' pdbref=',pdbref
 311
 312 C Read sequence if not taken from the pdb file.
 313       if (iscode.gt.0) then
 314         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 315       else
 316         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 317       endif
 318 C Convert sequence to numeric code
 319       do i=1,nres
 320         itype(i)=rescode(i,sequence(i),iscode)
 321       enddo
 322       print *,nres
 323       print '(20i4)',(itype(i),i=1,nres)
 324
 325       do i=1,nres
 326 #ifdef PROCOR
 327         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 328 #else
 329         if (itype(i).eq.ntyp1) then
 330 #endif
 331           itel(i)=0
 332 #ifdef PROCOR
 333         else if (iabs(itype(i+1)).ne.20) then
 334 #else
 335         else if (iabs(itype(i)).ne.20) then
 336 #endif
 337           itel(i)=1
 338         else
 339           itel(i)=2
 340         endif
 341       enddo
 342       write (iout,*) "ITEL"
 343       do i=1,nres-1
 344         write (iout,*) i,itype(i),itel(i)
 345       enddo
 346
 347       print *,'Call Read_Bridge.'
 348       call read_bridge
 349 C this fragment reads diheadral constrains
 350       if (with_dihed_constr) then
 351
 352       read (inp,*) ndih_constr
 353       if (ndih_constr.gt.0) then
 354 C        read (inp,*) ftors
 355 C        write (iout,*) 'FTORS',ftors
 356 C ftors is the force constant for torsional quartic constrains
 357         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 358      &                i=1,ndih_constr)
 359         write (iout,*)
 360      &   'There are',ndih_constr,' constraints on phi angles.'
 361         do i=1,ndih_constr
 362           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 363      &  ftors(i)
 364         enddo
 365         do i=1,ndih_constr
 366           phi0(i)=deg2rad*phi0(i)
 367           drange(i)=deg2rad*drange(i)
 368         enddo
 369       endif ! endif ndif_constr.gt.0
 370       endif ! with_dihed_constr
 371       if (with_theta_constr) then
 372 C with_theta_constr is keyword allowing for occurance of theta constrains
 373       read (inp,*) ntheta_constr
 374 C ntheta_constr is the number of theta constrains
 375       if (ntheta_constr.gt.0) then
 376 C        read (inp,*) ftors
 377         read (inp,*) (itheta_constr(i),theta_constr0(i),
 378      &  theta_drange(i),for_thet_constr(i),
 379      &  i=1,ntheta_constr)
 380 C the above code reads from 1 to ntheta_constr
 381 C itheta_constr(i) residue i for which is theta_constr
 382 C theta_constr0 the global minimum value
 383 C theta_drange is range for which there is no energy penalty
 384 C for_thet_constr is the force constant for quartic energy penalty
 385 C E=k*x**4
 386 C        if(me.eq.king.or..not.out1file)then
 387          write (iout,*)
 388      &   'There are',ntheta_constr,' constraints on phi angles.'
 389          do i=1,ntheta_constr
 390           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 391      &    theta_drange(i),
 392      &    for_thet_constr(i)
 393          enddo
 394 C        endif
 395         do i=1,ntheta_constr
 396           theta_constr0(i)=deg2rad*theta_constr0(i)
 397           theta_drange(i)=deg2rad*theta_drange(i)
 398         enddo
 399 C        if(me.eq.king.or..not.out1file)
 400 C     &   write (iout,*) 'FTORS',ftors
 401 C        do i=1,ntheta_constr
 402 C          ii = itheta_constr(i)
 403 C          thetabound(1,ii) = phi0(i)-drange(i)
 404 C          thetabound(2,ii) = phi0(i)+drange(i)
 405 C        enddo
 406       endif ! ntheta_constr.gt.0
 407       endif! with_theta_constr
 408
 409       nnt=1
 410       nct=nres
 411       print *,'NNT=',NNT,' NCT=',NCT
 412       if (itype(1).eq.ntyp1) nnt=2
 413       if (itype(nres).eq.ntyp1) nct=nct-1
 414       if (nstart.lt.nnt) nstart=nnt
 415       if (nend.gt.nct .or. nend.eq.0) nend=nct
 416       write (iout,*) "nstart",nstart," nend",nend
 417       nres0=nres
 418       if (constr_homology.gt.0) then
 419         call read_constr_homology(print_homology_restraints)
 420       endif
 421
 422 c      if (pdbref) then
 423 c        read(inp,'(a)') pdbfile
 424 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 425 c        open(ipdbin,file=pdbfile,status='old',err=33)
 426 c        goto 34
 427 c  33    write (iout,'(a)') 'Error opening PDB file.'
 428 c        stop
 429 c  34    continue
 430 c        print *,'Begin reading pdb data'
 431 c        call readpdb
 432 c        print *,'Finished reading pdb data'
 433 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 434 c        do i=1,nres
 435 c          itype_pdb(i)=itype(i)
 436 c        enddo
 437 c        close (ipdbin)
 438 c        write (iout,'(a,i3)') 'nsup=',nsup
 439 c        nstart_seq=nnt
 440 c        if (nsup.le.(nct-nnt+1)) then
 441 c          do i=0,nct-nnt+1-nsup
 442 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 443 c              nstart_seq=nnt+i
 444 c              goto 111
 445 c            endif
 446 c          enddo
 447 c          write (iout,'(a)')
 448 c     &            'Error - sequences to be superposed do not match.'
 449 c          stop
 450 c        else
 451 c          do i=0,nsup-(nct-nnt+1)
 452 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 453 c     &      then
 454 c              nstart_sup=nstart_sup+i
 455 c              nsup=nct-nnt+1
 456 c              goto 111
 457 c            endif
 458 c          enddo
 459 c          write (iout,'(a)')
 460 c     &            'Error - sequences to be superposed do not match.'
 461 c        endif
 462 c  111   continue
 463 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 464 c     &                 ' nstart_seq=',nstart_seq
 465 c      endif
 466       call init_int_table
 467       call setup_var
 468       write (iout,*) "molread: REFSTR",refstr
 469       if (refstr) then
 470         if (.not.pdbref) then
 471           call read_angles(inp,*38)
 472           goto 39
 473    38     write (iout,'(a)') 'Error reading reference structure.'
 474 #ifdef MPL
 475           call mp_stopall(Error_Msg)
 476 #else
 477           stop 'Error reading reference structure'
 478 #endif
 479    39     call chainbuild
 480           nstart_sup=nnt
 481           nstart_seq=nnt
 482           nsup=nct-nnt+1
 483           kkk=1
 484           do i=1,2*nres
 485             do j=1,3
 486               cref(j,i,kkk)=c(j,i)
 487             enddo
 488           enddo
 489         endif
 490         call contact(.true.,ncont_ref,icont_ref)
 491       endif
 492        if (ns.gt.0) then
 493 C        write (iout,'(/a,i3,a)')
 494 C     &  'The chain contains',ns,' disulfide-bridging cysteines.'
 495         write (iout,'(20i4)') (iss(i),i=1,ns)
 496        if (dyn_ss) then
 497           write(iout,*)"Running with dynamic disulfide-bond formation"
 498        else
 499         write (iout,'(/a/)') 'Pre-formed links are:'
 500         do i=1,nss
 501           i1=ihpb(i)-nres
 502           i2=jhpb(i)-nres
 503           it1=itype(i1)
 504           it2=itype(i2)
 505           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 506      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 507      &    ebr,forcon(i)
 508         enddo
 509         write (iout,'(a)')
 510        endif
 511       endif
 512       if (ns.gt.0.and.dyn_ss) then
 513           do i=nss+1,nhpb
 514             ihpb(i-nss)=ihpb(i)
 515             jhpb(i-nss)=jhpb(i)
 516             forcon(i-nss)=forcon(i)
 517             dhpb(i-nss)=dhpb(i)
 518           enddo
 519           nhpb=nhpb-nss
 520           nss=0
 521           call hpb_partition
 522           do i=1,ns
 523             dyn_ss_mask(iss(i))=.true.
 524           enddo
 525       endif
 526 c Read distance restraints
 527       if (constr_dist.gt.0) then
 528         call read_dist_constr
 529         call hpb_partition
 530       endif
 531       return
 532       end
 533 c-----------------------------------------------------------------------------
 534       logical function seq_comp(itypea,itypeb,length)
 535       implicit none
 536       integer length,itypea(length),itypeb(length)
 537       integer i
 538       do i=1,length
 539         if (itypea(i).ne.itypeb(i)) then
 540           seq_comp=.false.
 541           return
 542         endif
 543       enddo
 544       seq_comp=.true.
 545       return
 546       end
 547 c-----------------------------------------------------------------------------
 548       subroutine read_bridge
 549 C Read information about disulfide bridges.
 550       implicit none
 551       include 'DIMENSIONS'
 552       include 'COMMON.IOUNITS'
 553       include 'COMMON.GEO'
 554       include 'COMMON.VAR'
 555       include 'COMMON.INTERACT'
 556       include 'COMMON.LOCAL'
 557       include 'COMMON.NAMES'
 558       include 'COMMON.CHAIN'
 559       include 'COMMON.FFIELD'
 560       include 'COMMON.SBRIDGE'
 561       include 'COMMON.HEADER'
 562       include 'COMMON.CONTROL'
 563       include 'COMMON.TIME1'
 564 #ifdef MPL
 565       include 'COMMON.INFO'
 566 #endif
 567       integer i,j
 568 C Read bridging residues.
 569       read (inp,*) ns,(iss(i),i=1,ns)
 570       print *,'ns=',ns
 571 C Check whether the specified bridging residues are cystines.
 572       do i=1,ns
 573         if (itype(iss(i)).ne.1) then
 574           write (iout,'(2a,i3,a)')
 575      &   'Do you REALLY think that the residue ',
 576      &    restyp(itype(iss(i))),i,
 577      &   ' can form a disulfide bridge?!!!'
 578           write (*,'(2a,i3,a)')
 579      &   'Do you REALLY think that the residue ',
 580      &   restyp(itype(iss(i))),i,
 581      &   ' can form a disulfide bridge?!!!'
 582 #ifdef MPL
 583          call mp_stopall(error_msg)
 584 #else
 585          stop
 586 #endif
 587         endif
 588       enddo
 589 C Read preformed bridges.
 590       if (ns.gt.0) then
 591       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 592       if (nss.gt.0) then
 593         nhpb=nss
 594 C Check if the residues involved in bridges are in the specified list of
 595 C bridging residues.
 596         do i=1,nss
 597           do j=1,i-1
 598             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 599      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 600               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 601      &      ' contains residues present in other pairs.'
 602               write (*,'(a,i3,a)') 'Disulfide pair',i,
 603      &      ' contains residues present in other pairs.'
 604 #ifdef MPL
 605               call mp_stopall(error_msg)
 606 #else
 607               stop
 608 #endif
 609             endif
 610           enddo
 611           do j=1,ns
 612             if (ihpb(i).eq.iss(j)) goto 10
 613           enddo
 614           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 615    10     continue
 616           do j=1,ns
 617             if (jhpb(i).eq.iss(j)) goto 20
 618           enddo
 619           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 620    20     continue
 621 C          dhpb(i)=dbr
 622 C          forcon(i)=fbr
 623         enddo
 624         do i=1,nss
 625           ihpb(i)=ihpb(i)+nres
 626           jhpb(i)=jhpb(i)+nres
 627         enddo
 628       endif
 629       endif
 630       return
 631       end
 632 c----------------------------------------------------------------------------
 633       subroutine read_angles(kanal,*)
 634       implicit none
 635       include 'DIMENSIONS'
 636       include 'COMMON.GEO'
 637       include 'COMMON.VAR'
 638       include 'COMMON.CHAIN'
 639       include 'COMMON.IOUNITS'
 640       integer i,kanal
 641       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 642       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 643       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 644       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 645       do i=1,nres
 646         theta(i)=deg2rad*theta(i)
 647         phi(i)=deg2rad*phi(i)
 648         alph(i)=deg2rad*alph(i)
 649         omeg(i)=deg2rad*omeg(i)
 650       enddo
 651       return
 652    10 return1
 653       end
 654 c----------------------------------------------------------------------------
 655       subroutine reada(rekord,lancuch,wartosc,default)
 656       implicit none
 657       character*(*) rekord,lancuch
 658       double precision wartosc,default
 659       integer ilen,iread
 660       external ilen
 661       iread=index(rekord,lancuch)
 662       if (iread.eq.0) then
 663         wartosc=default
 664         return
 665       endif
 666       iread=iread+ilen(lancuch)+1
 667       read (rekord(iread:),*) wartosc
 668       return
 669       end
 670 c----------------------------------------------------------------------------
 671       subroutine multreada(rekord,lancuch,tablica,dim,default)
 672       implicit none
 673       integer dim,i
 674       double precision tablica(dim),default
 675       character*(*) rekord,lancuch
 676       integer ilen,iread
 677       external ilen
 678       do i=1,dim
 679         tablica(i)=default
 680       enddo
 681       iread=index(rekord,lancuch)
 682       if (iread.eq.0) return
 683       iread=iread+ilen(lancuch)+1
 684       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 685    10 return
 686       end
 687 c----------------------------------------------------------------------------
 688       subroutine readi(rekord,lancuch,wartosc,default)
 689       implicit none
 690       character*(*) rekord,lancuch
 691       integer wartosc,default
 692       integer ilen,iread
 693       external ilen
 694       iread=index(rekord,lancuch)
 695       if (iread.eq.0) then
 696         wartosc=default
 697         return
 698       endif
 699       iread=iread+ilen(lancuch)+1
 700       read (rekord(iread:),*) wartosc
 701       return
 702       end
 703 C----------------------------------------------------------------------
 704       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 705       implicit none
 706       integer dim,i
 707       integer tablica(dim),default
 708       character*(*) rekord,lancuch
 709       character*80 aux
 710       integer ilen,iread
 711       external ilen
 712       do i=1,dim
 713         tablica(i)=default
 714       enddo
 715       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 716       if (iread.eq.0) return
 717       iread=iread+ilen(lancuch)+1
 718       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 719    10 return
 720       end
 721
 722 c----------------------------------------------------------------------------
 723       subroutine card_concat(card)
 724       include 'DIMENSIONS'
 725       include 'COMMON.IOUNITS'
 726       character*(*) card
 727       character*80 karta,ucase
 728       external ilen
 729       read (inp,'(a)') karta
 730       karta=ucase(karta)
 731       card=' '
 732       do while (karta(80:80).eq.'&')
 733         card=card(:ilen(card)+1)//karta(:79)
 734         read (inp,'(a)') karta
 735         karta=ucase(karta)
 736       enddo
 737       card=card(:ilen(card)+1)//karta
 738       return
 739       end
 740 c----------------------------------------------------------------------------
 741       subroutine openunits
 742       implicit none
 743       include 'DIMENSIONS'
 744 #ifdef MPI
 745       include "mpif.h"
 746       character*3 liczba
 747       include "COMMON.MPI"
 748 #endif
 749       include 'COMMON.IOUNITS'
 750       include 'COMMON.CONTROL'
 751       integer lenpre,lenpot,ilen
 752       external ilen
 753       character*16 cformat,cprint
 754       character*16 ucase
 755       integer lenint,lenout
 756       call getenv('INPUT',prefix)
 757       call getenv('OUTPUT',prefout)
 758       call getenv('INTIN',prefintin)
 759       call getenv('COORD',cformat)
 760       call getenv('PRINTCOOR',cprint)
 761       call getenv('SCRATCHDIR',scratchdir)
 762       from_bx=.true.
 763       from_cx=.false.
 764       if (index(ucase(cformat),'CX').gt.0) then
 765         from_cx=.true.
 766         from_bx=.false.
 767       endif
 768       from_cart=.true.
 769       lenpre=ilen(prefix)
 770       lenout=ilen(prefout)
 771       lenint=ilen(prefintin)
 772 C Get the names and open the input files
 773       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 774 #ifdef MPI
 775       write (liczba,'(bz,i3.3)') me
 776       outname=prefout(:lenout)//'_clust.out_'//liczba
 777 #else
 778       outname=prefout(:lenout)//'_clust.out'
 779 #endif
 780       if (from_bx) then
 781         intinname=prefintin(:lenint)//'.bx'
 782       else if (from_cx) then
 783         intinname=prefintin(:lenint)//'.cx'
 784       else
 785         intinname=prefintin(:lenint)//'.int'
 786       endif
 787       rmsname=prefintin(:lenint)//'.rms'
 788       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 789      &       status='unknown')
 790       open (jrms,file=rmsname,status='unknown')
 791       open(iout,file=outname,status='unknown')
 792 C Get parameter filenames and open the parameter files.
 793       call getenv('BONDPAR',bondname)
 794       open (ibond,file=bondname,status='old')
 795       call getenv('THETPAR',thetname)
 796       open (ithep,file=thetname,status='old')
 797       call getenv('ROTPAR',rotname)
 798       open (irotam,file=rotname,status='old')
 799       call getenv('TORPAR',torname)
 800       open (itorp,file=torname,status='old')
 801       call getenv('TORDPAR',tordname)
 802       open (itordp,file=tordname,status='old')
 803       call getenv('FOURIER',fouriername)
 804       open (ifourier,file=fouriername,status='old')
 805       call getenv('ELEPAR',elename)
 806       open (ielep,file=elename,status='old')
 807       call getenv('SIDEPAR',sidename)
 808       open (isidep,file=sidename,status='old')
 809       call getenv('SIDEP',sidepname)
 810       open (isidep1,file=sidepname,status="old")
 811       call getenv('SCCORPAR',sccorname)
 812       open (isccor,file=sccorname,status="old")
 813       call getenv('LIPTRANPAR',liptranname)
 814       open (iliptranpar,file=liptranname,status='old')
 815 #ifndef OLDSCP
 816 C
 817 C 8/9/01 In the newest version SCp interaction constants are read from a file
 818 C Use -DOLDSCP to use hard-coded constants instead.
 819 C
 820       call getenv('SCPPAR',scpname)
 821       open (iscpp,file=scpname,status='old')
 822 #endif
 823       return
 824       end
 825       subroutine read_dist_constr
 826       implicit real*8 (a-h,o-z)
 827       include 'DIMENSIONS'
 828 #ifdef MPI
 829       include 'mpif.h'
 830 #endif
 831       include 'COMMON.CONTROL'
 832       include 'COMMON.CHAIN'
 833       include 'COMMON.IOUNITS'
 834       include 'COMMON.SBRIDGE'
 835       integer ifrag_(2,100),ipair_(2,100)
 836       double precision wfrag_(100),wpair_(100)
 837       character*500 controlcard
 838       logical lprn /.true./
 839       write (iout,*) "Calling read_dist_constr"
 840 C      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 841 C      call flush(iout)
 842       write(iout,*) "TU sie wywalam?"
 843       call card_concat(controlcard)
 844       write (iout,*) controlcard
 845       call flush(iout)
 846       call readi(controlcard,"NFRAG",nfrag_,0)
 847       call readi(controlcard,"NPAIR",npair_,0)
 848       call readi(controlcard,"NDIST",ndist_,0)
 849       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 850       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 851       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 852       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 853       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 854       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 855       write (iout,*) "IFRAG"
 856       do i=1,nfrag_
 857         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 858       enddo
 859       write (iout,*) "IPAIR"
 860       do i=1,npair_
 861         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 862       enddo
 863       call flush(iout)
 864       do i=1,nfrag_
 865         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 866         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 867      &    ifrag_(2,i)=nstart_sup+nsup-1
 868 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 869         call flush(iout)
 870         if (wfrag_(i).gt.0.0d0) then
 871         do j=ifrag_(1,i),ifrag_(2,i)-1
 872           do k=j+1,ifrag_(2,i)
 873             write (iout,*) "j",j," k",k
 874             ddjk=dist(j,k)
 875             if (constr_dist.eq.1) then
 876               nhpb=nhpb+1
 877               ihpb(nhpb)=j
 878               jhpb(nhpb)=k
 879               dhpb(nhpb)=ddjk
 880               forcon(nhpb)=wfrag_(i)
 881             else if (constr_dist.eq.2) then
 882               if (ddjk.le.dist_cut) then
 883                 nhpb=nhpb+1
 884                 ihpb(nhpb)=j
 885                 jhpb(nhpb)=k
 886                 dhpb(nhpb)=ddjk
 887                 forcon(nhpb)=wfrag_(i)
 888               endif
 889             else
 890               nhpb=nhpb+1
 891               ihpb(nhpb)=j
 892               jhpb(nhpb)=k
 893               dhpb(nhpb)=ddjk
 894               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 895             endif
 896             if (lprn)
 897      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 898      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 899           enddo
 900         enddo
 901         endif
 902       enddo
 903       do i=1,npair_
 904         if (wpair_(i).gt.0.0d0) then
 905         ii = ipair_(1,i)
 906         jj = ipair_(2,i)
 907         if (ii.gt.jj) then
 908           itemp=ii
 909           ii=jj
 910           jj=itemp
 911         endif
 912         do j=ifrag_(1,ii),ifrag_(2,ii)
 913           do k=ifrag_(1,jj),ifrag_(2,jj)
 914             nhpb=nhpb+1
 915             ihpb(nhpb)=j
 916             jhpb(nhpb)=k
 917             forcon(nhpb)=wpair_(i)
 918             dhpb(nhpb)=dist(j,k)
 919             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 920      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 921           enddo
 922         enddo
 923         endif
 924       enddo
 925       do i=1,ndist_
 926         if (constr_dist.eq.11) then
 927         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 928      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 929         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 930 C        write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 931 C     &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 932         else
 933         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 934         endif
 935         if (forcon(nhpb+1).gt.0.0d0) then
 936           nhpb=nhpb+1
 937           if (ibecarb(i).gt.0) then
 938             ihpb(i)=ihpb(i)+nres
 939             jhpb(i)=jhpb(i)+nres
 940           endif
 941           if (dhpb(nhpb).eq.0.0d0)
 942      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 943 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 944           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 945      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 946         endif
 947 C      endif
 948       enddo
 949       call hpb_partition
 950 #ifdef AIX
 951       call flush_(iout)
 952 #else
 953       call flush(iout)
 954 #endif
 955       return
 956       end
 957
 958 c====-------------------------------------------------------------------
 959       subroutine read_constr_homology(lprn)
 960
 961       include 'DIMENSIONS'
 962 #ifdef MPI
 963       include 'mpif.h'
 964 #endif
 965       include 'COMMON.SETUP'
 966       include 'COMMON.CONTROL'
 967       include 'COMMON.CHAIN'
 968       include 'COMMON.IOUNITS'
 969       include 'COMMON.GEO'
 970       include 'COMMON.INTERACT'
 971       include 'COMMON.NAMES'
 972       include 'COMMON.HOMRESTR'
 973 c
 974 c For new homol impl
 975 c
 976       include 'COMMON.VAR'
 977 c     include 'include_unres/COMMON.VAR'
 978 c
 979
 980 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 981 c    &                 dist_cut
 982 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 983 c    &    sigma_odl_temp(maxres,maxres,max_template)
 984       character*2 kic2
 985       character*24 model_ki_dist, model_ki_angle
 986       character*500 controlcard
 987       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 988       integer idomain(max_template,maxres)
 989       integer ilen
 990       external ilen
 991       logical lprn
 992       logical unres_pdb,liiflag
 993 c
 994 c     FP - Nov. 2014 Temporary specifications for new vars
 995 c
 996       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
 997      &                 rescore3_tmp
 998       double precision, dimension (max_template,maxres) :: rescore
 999       double precision, dimension (max_template,maxres) :: rescore2
1000       double precision, dimension (max_template,maxres) :: rescore3
1001       character*24 tpl_k_rescore
1002 c -----------------------------------------------------------------
1003 c Reading multiple PDB ref structures and calculation of retraints
1004 c not using pre-computed ones stored in files model_ki_{dist,angle}
1005 c FP (Nov., 2014)
1006 c -----------------------------------------------------------------
1007 c
1008 c
1009 c Alternative: reading from input
1010 #ifdef DEBUG
1011       write (iout,*) "BEGIN READ HOMOLOGY INFO"
1012 #ifdef AIX
1013       call flush_(iout)
1014 #else
1015       call flush(iout)
1016 #endif
1017 #endif
1018       call card_concat(controlcard)
1019       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1020       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1021       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1022       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1023       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1024       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1025       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1026       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1027       if (homol_nset.gt.1)then
1028          call readi(controlcard,"ISET",iset,1)
1029          call card_concat(controlcard)
1030          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1031       else
1032         iset=1
1033         waga_homology(1)=1.0
1034       endif
1035 c
1036 #ifdef DEBUG
1037       write(iout,*) "read_constr_homology iset",iset
1038       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1039 #ifdef AIX
1040       call flush_(iout)
1041 #else
1042       call flush(iout)
1043 #endif
1044 #endif
1045 cd      write (iout,*) "nnt",nnt," nct",nct
1046 cd      call flush(iout)
1047
1048
1049       lim_odl=0
1050       lim_dih=0
1051 c
1052 c  New
1053 c
1054 c
1055 c  Reading HM global scores (prob not required)
1056 c
1057       do i = nnt,nct
1058         do k=1,constr_homology
1059          idomain(k,i)=0
1060         enddo
1061       enddo
1062 c     open (4,file="HMscore")
1063 c     do k=1,constr_homology
1064 c       read (4,*,end=521) hmscore_tmp
1065 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1066 c       write(*,*) "Model", k, ":", hmscore(k)
1067 c     enddo
1068 c521  continue
1069
1070       ii=0
1071       do i = nnt,nct-2
1072         do j=i+2,nct
1073         ii=ii+1
1074         ii_in_use(ii)=0
1075         enddo
1076       enddo
1077 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1078
1079       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1080       do k=1,constr_homology
1081
1082         read(inp,'(a)') pdbfile
1083 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1084 c  Next stament causes error upon compilation (?)
1085 c       if(me.eq.king.or. .not. out1file)
1086 c         write (iout,'(2a)') 'PDB data will be read from file ',
1087 c    &   pdbfile(:ilen(pdbfile))
1088          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1089      &  pdbfile(:ilen(pdbfile))
1090         open(ipdbin,file=pdbfile,status='old',err=33)
1091         goto 34
1092   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1093      &  pdbfile(:ilen(pdbfile))
1094         stop
1095   34    continue
1096 c        print *,'Begin reading pdb data'
1097 c
1098 c Files containing res sim or local scores (former containing sigmas)
1099 c
1100
1101         write(kic2,'(bz,i2.2)') k
1102
1103         tpl_k_rescore="template"//kic2//".sco"
1104
1105         unres_pdb=.false.
1106         call readpdb_template(k)
1107         do i=1,2*nres
1108           do j=1,3
1109             crefjlee(j,i)=c(j,i)
1110           enddo
1111         enddo
1112 #ifdef DEBUG
1113         do i=1,nres
1114           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1115      &      (crefjlee(j,i+nres),j=1,3)
1116         enddo
1117         write (iout,*) "READ HOMOLOGY INFO"
1118         write (iout,*) "read_constr_homology x: after reading pdb file"
1119         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1120         write (iout,*) "waga_dist",waga_dist
1121         write (iout,*) "waga_angle",waga_angle
1122         write (iout,*) "waga_theta",waga_theta
1123         write (iout,*) "waga_d",waga_d
1124         write (iout,*) "dist_cut",dist_cut
1125 #endif
1126 #ifdef AIX
1127       call flush_(iout)
1128 #else
1129       call flush(iout)
1130 #endif
1131
1132 c
1133 c     Distance restraints
1134 c
1135 c          ... --> odl(k,ii)
1136 C Copy the coordinates from reference coordinates (?)
1137         do i=1,2*nres
1138           do j=1,3
1139 c            c(j,i)=cref(j,i)
1140 c           write (iout,*) "c(",j,i,") =",c(j,i)
1141           enddo
1142         enddo
1143 c
1144 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1145 c
1146 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1147           open (ientin,file=tpl_k_rescore,status='old')
1148           if (nnt.gt.1) rescore(k,1)=0.0d0
1149           do irec=nnt,nct ! loop for reading res sim
1150             if (read2sigma) then
1151              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1152      &                                rescore3_tmp,idomain_tmp
1153              i_tmp=i_tmp+nnt-1
1154              idomain(k,i_tmp)=idomain_tmp
1155              rescore(k,i_tmp)=rescore_tmp
1156              rescore2(k,i_tmp)=rescore2_tmp
1157              rescore3(k,i_tmp)=rescore3_tmp
1158              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1159      &                      i_tmp,rescore2_tmp,rescore_tmp,
1160      &                                rescore3_tmp,idomain_tmp
1161             else
1162              idomain(k,irec)=1
1163              read (ientin,*,end=1401) rescore_tmp
1164
1165 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1166              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1167 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1168             endif
1169           enddo
1170  1401   continue
1171         close (ientin)
1172         if (waga_dist.ne.0.0d0) then
1173           ii=0
1174           do i = nnt,nct-2
1175             do j=i+2,nct
1176
1177               x12=c(1,i)-c(1,j)
1178               y12=c(2,i)-c(2,j)
1179               z12=c(3,i)-c(3,j)
1180               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1181 c              write (iout,*) k,i,j,distal,dist2_cut
1182
1183             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1184      &            .and. distal.le.dist2_cut ) then
1185
1186               ii=ii+1
1187               ii_in_use(ii)=1
1188               l_homo(k,ii)=.true.
1189
1190 c             write (iout,*) "k",k
1191 c             write (iout,*) "i",i," j",j," constr_homology",
1192 c    &                       constr_homology
1193               ires_homo(ii)=i
1194               jres_homo(ii)=j
1195               odl(k,ii)=distal
1196               if (read2sigma) then
1197                 sigma_odl(k,ii)=0
1198                 do ik=i,j
1199                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1200                 enddo
1201                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1202                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1203      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1204               else
1205                 if (odl(k,ii).le.dist_cut) then
1206                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1207                 else
1208 #ifdef OLDSIGMA
1209                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1210      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1211 #else
1212                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1213      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1214 #endif
1215                 endif
1216               endif
1217               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1218             else
1219               ii=ii+1
1220               l_homo(k,ii)=.false.
1221             endif
1222             enddo
1223           enddo
1224         lim_odl=ii
1225         endif
1226 c
1227 c     Theta, dihedral and SC retraints
1228 c
1229         if (waga_angle.gt.0.0d0) then
1230 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1231 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1232 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1233 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1234 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1235 c    &                            sigma_dih(k,i+nnt-1)
1236 c         enddo
1237 c1402   continue
1238 c         close (ientin)
1239           do i = nnt+3,nct
1240             if (idomain(k,i).eq.0) then
1241                sigma_dih(k,i)=0.0
1242                cycle
1243             endif
1244             dih(k,i)=phiref(i) ! right?
1245 c           read (ientin,*) sigma_dih(k,i) ! original variant
1246 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1247 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1248 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1249 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1250 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1251
1252             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1253      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1254 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1255 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1256 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1257 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1258 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1259 c   Instead of res sim other local measure of b/b str reliability possible
1260             if (sigma_dih(k,i).ne.0)
1261      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1262 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1263           enddo
1264           lim_dih=nct-nnt-2
1265         endif
1266
1267         if (waga_theta.gt.0.0d0) then
1268 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1269 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1270 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1271 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1272 c    &                              sigma_theta(k,i+nnt-1)
1273 c         enddo
1274 c1403   continue
1275 c         close (ientin)
1276
1277           do i = nnt+2,nct ! right? without parallel.
1278 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1279 c         do i=ithet_start,ithet_end ! with FG parallel.
1280              if (idomain(k,i).eq.0) then
1281               sigma_theta(k,i)=0.0
1282               cycle
1283              endif
1284              thetatpl(k,i)=thetaref(i)
1285 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1286 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1287 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1288 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1289 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1290              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1291      &                        rescore(k,i-2))/3.0
1292 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1293              if (sigma_theta(k,i).ne.0)
1294      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1295
1296 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1297 c                             rescore(k,i-2) !  right expression ?
1298 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1299           enddo
1300         endif
1301
1302         if (waga_d.gt.0.0d0) then
1303 c       open (ientin,file=tpl_k_sigma_d,status='old')
1304 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1305 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1306 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1307 c    &                          sigma_d(k,i+nnt-1)
1308 c         enddo
1309 c1404   continue
1310
1311           do i = nnt,nct ! right? without parallel.
1312 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1313 c         do i=loc_start,loc_end ! with FG parallel.
1314                if (itype(i).eq.10) cycle
1315                if (idomain(k,i).eq.0 ) then
1316                   sigma_d(k,i)=0.0
1317                   cycle
1318                endif
1319                xxtpl(k,i)=xxref(i)
1320                yytpl(k,i)=yyref(i)
1321                zztpl(k,i)=zzref(i)
1322 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1323 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1324 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1325 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1326                sigma_d(k,i)=rescore3(k,i) !  right expression ?
1327                if (sigma_d(k,i).ne.0)
1328      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1329
1330 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1331 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1332 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1333           enddo
1334         endif
1335       enddo
1336 c
1337 c remove distance restraints not used in any model from the list
1338 c shift data in all arrays
1339 c
1340       if (waga_dist.ne.0.0d0) then
1341         ii=0
1342         liiflag=.true.
1343         do i=nnt,nct-2
1344          do j=i+2,nct
1345           ii=ii+1
1346           if (ii_in_use(ii).eq.0.and.liiflag) then
1347             liiflag=.false.
1348             iistart=ii
1349           endif
1350           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1351      &                   .not.liiflag.and.ii.eq.lim_odl) then
1352              if (ii.eq.lim_odl) then
1353               iishift=ii-iistart+1
1354              else
1355               iishift=ii-iistart
1356              endif
1357              liiflag=.true.
1358              do ki=iistart,lim_odl-iishift
1359               ires_homo(ki)=ires_homo(ki+iishift)
1360               jres_homo(ki)=jres_homo(ki+iishift)
1361               ii_in_use(ki)=ii_in_use(ki+iishift)
1362               do k=1,constr_homology
1363                odl(k,ki)=odl(k,ki+iishift)
1364                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1365                l_homo(k,ki)=l_homo(k,ki+iishift)
1366               enddo
1367              enddo
1368              ii=ii-iishift
1369              lim_odl=lim_odl-iishift
1370           endif
1371          enddo
1372         enddo
1373       endif
1374       if (constr_homology.gt.0) call homology_partition
1375       if (constr_homology.gt.0) call init_int_table
1376 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1377 cd     & "lim_xx=",lim_xx
1378 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1379 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1380 c
1381 c Print restraints
1382 c
1383       if (.not.lprn) return
1384 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1385       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1386        write (iout,*) "Distance restraints from templates"
1387        do ii=1,lim_odl
1388        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1389      &  ii,ires_homo(ii),jres_homo(ii),
1390      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1391      &  ki=1,constr_homology)
1392        enddo
1393        write (iout,*) "Dihedral angle restraints from templates"
1394        do i=nnt+3,nct
1395         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1396      &      (rad2deg*dih(ki,i),
1397      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1398        enddo
1399        write (iout,*) "Virtual-bond angle restraints from templates"
1400        do i=nnt+2,nct
1401         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1402      &      (rad2deg*thetatpl(ki,i),
1403      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1404        enddo
1405        write (iout,*) "SC restraints from templates"
1406        do i=nnt,nct
1407         write(iout,'(i5,100(4f8.2,4x))') i,
1408      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1409      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1410        enddo
1411       endif
1412 c -----------------------------------------------------------------
1413       return
1414       end