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