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
 936       logical lprn
 937       logical unres_pdb
 938 c
 939 c     FP - Nov. 2014 Temporary specifications for new vars
 940 c
 941       double precision rescore_tmp,x12,y12,z12
 942       double precision, dimension (max_template,maxres) :: rescore
 943       character*24 tpl_k_rescore
 944 c -----------------------------------------------------------------
 945 c Reading multiple PDB ref structures and calculation of retraints
 946 c not using pre-computed ones stored in files model_ki_{dist,angle}
 947 c FP (Nov., 2014)
 948 c -----------------------------------------------------------------
 949 c
 950 c
 951 c Alternative: reading from input
 952 #ifdef DEBUG
 953       write (iout,*) "BEGIN READ HOMOLOGY INFO"
 954 #ifdef AIX
 955       call flush_(iout)
 956 #else
 957       call flush(iout)
 958 #endif
 959 #endif
 960       call card_concat(controlcard)
 961       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 962       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 963       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 964       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 965       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 966
 967       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
 968       if (homol_nset.gt.1)then
 969          call readi(controlcard,"ISET",iset,1)
 970          call card_concat(controlcard)
 971          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 972       else
 973         iset=1
 974         waga_homology(1)=1.0
 975       endif
 976 c
 977 #ifdef DEBUG
 978       write(iout,*) "read_constr_homology iset",iset
 979       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 980 #ifdef AIX
 981       call flush_(iout)
 982 #else
 983       call flush(iout)
 984 #endif
 985 #endif
 986 cd      write (iout,*) "nnt",nnt," nct",nct
 987 cd      call flush(iout)
 988
 989
 990       lim_odl=0
 991       lim_dih=0
 992 c
 993 c  New
 994 c
 995       lim_theta=0
 996       lim_xx=0
 997 c
 998 c  Reading HM global scores (prob not required)
 999 c
1000 c     open (4,file="HMscore")
1001 c     do k=1,constr_homology
1002 c       read (4,*,end=521) hmscore_tmp
1003 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1004 c       write(*,*) "Model", k, ":", hmscore(k)
1005 c     enddo
1006 c521  continue
1007
1008 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1009
1010       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1011       do k=1,constr_homology
1012
1013         read(inp,'(a)') pdbfile
1014 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1015 c  Next stament causes error upon compilation (?)
1016 c       if(me.eq.king.or. .not. out1file)
1017 c         write (iout,'(2a)') 'PDB data will be read from file ',
1018 c    &   pdbfile(:ilen(pdbfile))
1019         open(ipdbin,file=pdbfile,status='old',err=33)
1020         goto 34
1021   33    write (iout,'(a)') 'Error opening PDB file.'
1022         stop
1023   34    continue
1024 c        print *,'Begin reading pdb data'
1025 c
1026 c Files containing res sim or local scores (former containing sigmas)
1027 c
1028
1029         write(kic2,'(bz,i2.2)') k
1030
1031         tpl_k_rescore="template"//kic2//".sco"
1032 c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
1033 c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
1034 c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
1035 c       tpl_k_sigma_d="template"//kic2//".sigma_d"
1036
1037         unres_pdb=.false.
1038         call readpdb
1039         do i=1,2*nres
1040           do j=1,3
1041             crefjlee(j,i)=c(j,i)
1042           enddo
1043         enddo
1044 #ifdef DEBUG
1045         do i=1,nres
1046           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1047      &      (crefjlee(j,i+nres),j=1,3)
1048         enddo
1049         write (iout,*) "READ HOMOLOGY INFO"
1050         write (iout,*) "read_constr_homology x: after reading pdb file"
1051         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1052         write (iout,*) "waga_dist",waga_dist
1053         write (iout,*) "waga_angle",waga_angle
1054         write (iout,*) "waga_theta",waga_theta
1055         write (iout,*) "waga_d",waga_d
1056         write (iout,*) "dist_cut",dist_cut
1057 #endif
1058 #ifdef AIX
1059       call flush_(iout)
1060 #else
1061       call flush(iout)
1062 #endif
1063
1064 c
1065 c     Distance restraints
1066 c
1067 c          ... --> odl(k,ii)
1068 C Copy the coordinates from reference coordinates (?)
1069         do i=1,2*nres
1070           do j=1,3
1071             c(j,i)=cref(j,i)
1072 c           write (iout,*) "c(",j,i,") =",c(j,i)
1073           enddo
1074         enddo
1075 c
1076 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1077 c
1078 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1079           open (ientin,file=tpl_k_rescore,status='old')
1080           do irec=1,maxdim ! loop for reading res sim
1081             if (irec.eq.1) then
1082                rescore(k,irec)=0.0d0
1083                goto 1301
1084             endif
1085             read (ientin,*,end=1401) rescore_tmp
1086 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1087             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1088 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1089  1301     continue
1090           enddo
1091  1401   continue
1092           close (ientin)
1093 c         open (ientin,file=tpl_k_sigma_odl,status='old')
1094 c         do irec=1,maxdim ! loop for reading sigma_odl
1095 c            read (ientin,*,end=1401) i, j,
1096 c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
1097 c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
1098 c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k)
1099 c         enddo
1100 c 1401   continue
1101 c         close (ientin)
1102         if (waga_dist.ne.0.0d0) then
1103           ii=0
1104           do i = nnt,nct-2 ! right? without parallel.
1105             do j=i+2,nct ! right?
1106 c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
1107 c           do j=i+2,nres ! ibid
1108 c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
1109 c           do j=i+2,nct ! ibid
1110               ii=ii+1
1111 c             write (iout,*) "k",k
1112 c             write (iout,*) "i",i," j",j," constr_homology",
1113 c    &                       constr_homology
1114               ires_homo(ii)=i
1115               jres_homo(ii)=j
1116 c
1117 c Attempt to replace dist(i,j) by its definition in ...
1118 c
1119               x12=c(1,i)-c(1,j)
1120               y12=c(2,i)-c(2,j)
1121               z12=c(3,i)-c(3,j)
1122               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1123               odl(k,ii)=distal
1124 c
1125 c             odl(k,ii)=dist(i,j)
1126 c             write (iout,*) "dist(",i,j,") =",dist(i,j)
1127 c             write (iout,*) "distal = ",distal
1128 c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
1129 c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1130 c    &                      "rescore(",k,j,") =",rescore(k,j)
1131 c
1132 c  Calculation of sigma from res sim
1133 c
1134 c             if (odl(k,ii).le.6.0d0) then
1135 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
1136 c    Other functional forms possible depending on odl(k,ii), eg.
1137 c
1138             if (odl(k,ii).le.dist_cut) then
1139               sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
1140 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
1141             else
1142 #ifdef OLDSIGMA
1143               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
1144      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1145 #else
1146               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
1147      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1148 #endif
1149 c   Following expr replaced by a positive exp argument
1150 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1151 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
1152
1153 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
1154 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
1155             endif
1156 c
1157               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
1158 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
1159 c
1160 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
1161 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
1162             enddo
1163 c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
1164           enddo
1165 c         lim_odl=ii
1166 c         if (constr_homology.gt.0) call homology_partition
1167         endif
1168 c
1169 c     Theta, dihedral and SC retraints
1170 c
1171         if (waga_angle.gt.0.0d0) then
1172 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1173 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1174 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1175 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1176 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1177 c    &                            sigma_dih(k,i+nnt-1)
1178 c         enddo
1179 c1402   continue
1180 c         close (ientin)
1181           do i = nnt+3,nct ! right? without parallel.
1182 c         do i=1,nres ! alternative for bounds acc to readpdb?
1183 c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
1184 c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
1185             dih(k,i)=phiref(i) ! right?
1186 c           read (ientin,*) sigma_dih(k,i) ! original variant
1187 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1188 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1189 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1190 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1191 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1192
1193             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1194      &                     rescore(k,i-2)+rescore(k,i-3))/4.0  !  right expression ?
1195 c
1196 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1197 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1198 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1199 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1200 c   Instead of res sim other local measure of b/b str reliability possible
1201             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1202 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1203             if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
1204 c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
1205           enddo
1206         endif
1207
1208         if (waga_theta.gt.0.0d0) then
1209 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1210 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1211 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1212 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1213 c    &                              sigma_theta(k,i+nnt-1)
1214 c         enddo
1215 c1403   continue
1216 c         close (ientin)
1217
1218           do i = nnt+2,nct ! right? without parallel.
1219 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1220 c         do i=ithet_start,ithet_end ! with FG parallel.
1221              thetatpl(k,i)=thetaref(i)
1222 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1223 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1224 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1225 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1226 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1227              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1228      &                        rescore(k,i-2))/3.0 !  right expression ?
1229              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1230
1231 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1232 c                             rescore(k,i-2) !  right expression ?
1233 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1234              if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
1235           enddo
1236         endif
1237
1238         if (waga_d.gt.0.0d0) then
1239 c       open (ientin,file=tpl_k_sigma_d,status='old')
1240 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1241 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1242 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1243 c    &                          sigma_d(k,i+nnt-1)
1244 c         enddo
1245 c1404   continue
1246           close (ientin)
1247
1248           do i = nnt,nct ! right? without parallel.
1249 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1250 c         do i=loc_start,loc_end ! with FG parallel.
1251              if (itype(i).eq.10) goto 1 ! right?
1252                xxtpl(k,i)=xxref(i)
1253                yytpl(k,i)=yyref(i)
1254                zztpl(k,i)=zzref(i)
1255 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1256 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1257 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1258 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1259                sigma_d(k,i)=rescore(k,i) !  right expression ?
1260                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1261
1262 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1263 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1264 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1265                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1266     1     continue
1267           enddo
1268         endif
1269         close(ientin)
1270       enddo
1271       if (waga_dist.ne.0.0d0) lim_odl=ii
1272       if (constr_homology.gt.0) call homology_partition
1273       if (constr_homology.gt.0) call init_int_table
1274 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1275 cd     & "lim_xx=",lim_xx
1276 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1277 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1278 c
1279 c Print restraints
1280 c
1281       if (.not.lprn) return
1282 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1283       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1284        write (iout,*) "Distance restraints from templates"
1285        do ii=1,lim_odl
1286        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
1287      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
1288        enddo
1289        write (iout,*) "Dihedral angle restraints from templates"
1290        do i=nnt+3,lim_dih
1291         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
1292      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1293        enddo
1294        write (iout,*) "Virtual-bond angle restraints from templates"
1295        do i=nnt+2,lim_theta
1296         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
1297      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1298        enddo
1299        write (iout,*) "SC restraints from templates"
1300        do i=nnt,lim_xx
1301         write(iout,'(i5,10(4f8.2,4x))') i,
1302      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1303      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1304        enddo
1305       endif
1306 c -----------------------------------------------------------------
1307       return
1308       end
1309
1310