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