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