source/cluster/wham/src-old/readrtns.F.org

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