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