source/wham/src-M/molread_zs.F

   1       subroutine molread(*)
   2 C
   3 C Read molecular data.
   4 C
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7       include 'DIMENSIONS.ZSCOPT'
   8       include 'DIMENSIONS.FREE'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.GEO'
  11       include 'COMMON.VAR'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.LOCAL'
  14       include 'COMMON.NAMES'
  15       include 'COMMON.CHAIN'
  16       include 'COMMON.FFIELD'
  17       include 'COMMON.SBRIDGE'
  18       include 'COMMON.TORCNSTR'
  19       include 'COMMON.CONTROL'
  20       character*4 sequence(maxres)
  21       integer rescode
  22       double precision x(maxvar)
  23       character*320 controlcard,ucase
  24       dimension itype_pdb(maxres)
  25       logical seq_comp
  26       call card_concat(controlcard,.true.)
  27       call reada(controlcard,'SCAL14',scal14,0.4d0)
  28       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
  29       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
  30       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
  31       r0_corr=cutoff_corr-delt_corr
  32       call readi(controlcard,"NRES",nres,0)
  33       iscode=index(controlcard,"ONE_LETTER")
  34       if (nres.le.0) then
  35         write (iout,*) "Error: no residues in molecule"
  36         return1
  37       endif
  38       if (nres.gt.maxres) then
  39         write (iout,*) "Error: too many residues",nres,maxres
  40       endif
  41       write(iout,*) 'nres=',nres
  42 C Read sequence of the protein
  43       if (iscode.gt.0) then
  44         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
  45       else
  46         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
  47       endif
  48 C Convert sequence to numeric code
  49       do i=1,nres
  50         itype(i)=rescode(i,sequence(i),iscode)
  51       enddo
  52       write (iout,*) "Numeric code:"
  53       write (iout,'(20i4)') (itype(i),i=1,nres)
  54       do i=1,nres-1
  55 #ifdef PROCOR
  56         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
  57 #else
  58         if (itype(i).eq.ntyp1) then
  59 #endif
  60           itel(i)=0
  61 #ifdef PROCOR
  62         else if (iabs(itype(i+1)).ne.20) then
  63 #else
  64         else if (iabs(itype(i)).ne.20) then
  65 #endif
  66           itel(i)=1
  67         else
  68           itel(i)=2
  69         endif
  70       enddo
  71        write (iout,*) "ITEL"
  72        do i=1,nres-1
  73          write (iout,*) i,itype(i),itel(i)
  74        enddo
  75       call read_bridge
  76
  77       if (with_dihed_constr) then
  78
  79       read (inp,*) ndih_constr
  80       if (ndih_constr.gt.0) then
  81 C        read (inp,*) ftors
  82 C        write (iout,*) 'FTORS',ftors
  83         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
  84      &   i=1,ndih_constr)
  85         write (iout,*)
  86      &   'There are',ndih_constr,' constraints on phi angles.'
  87         do i=1,ndih_constr
  88           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
  89      &    ftors(i)
  90         enddo
  91         do i=1,ndih_constr
  92           phi0(i)=deg2rad*phi0(i)
  93           drange(i)=deg2rad*drange(i)
  94         enddo
  95       endif
  96
  97       endif
  98       if (with_theta_constr) then
  99 C with_theta_constr is keyword allowing for occurance of theta constrains
 100       read (inp,*) ntheta_constr
 101 C ntheta_constr is the number of theta constrains
 102       if (ntheta_constr.gt.0) then
 103 C        read (inp,*) ftors
 104         read (inp,*) (itheta_constr(i),theta_constr0(i),
 105      &  theta_drange(i),for_thet_constr(i),
 106      &  i=1,ntheta_constr)
 107 C the above code reads from 1 to ntheta_constr
 108 C itheta_constr(i) residue i for which is theta_constr
 109 C theta_constr0 the global minimum value
 110 C theta_drange is range for which there is no energy penalty
 111 C for_thet_constr is the force constant for quartic energy penalty
 112 C E=k*x**4
 113 C        if(me.eq.king.or..not.out1file)then
 114          write (iout,*)
 115      &   'There are',ntheta_constr,' constraints on phi angles.'
 116          do i=1,ntheta_constr
 117           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 118      &    theta_drange(i),
 119      &    for_thet_constr(i)
 120          enddo
 121 C        endif
 122         do i=1,ntheta_constr
 123           theta_constr0(i)=deg2rad*theta_constr0(i)
 124           theta_drange(i)=deg2rad*theta_drange(i)
 125         enddo
 126 C        if(me.eq.king.or..not.out1file)
 127 C     &   write (iout,*) 'FTORS',ftors
 128 C        do i=1,ntheta_constr
 129 C          ii = itheta_constr(i)
 130 C          thetabound(1,ii) = phi0(i)-drange(i)
 131 C          thetabound(2,ii) = phi0(i)+drange(i)
 132 C        enddo
 133       endif ! ntheta_constr.gt.0
 134       endif! with_theta_constr
 135       nnt=1
 136       nct=nres
 137       if (itype(1).eq.ntyp1) nnt=2
 138       if (itype(nres).eq.ntyp1) nct=nct-1
 139       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 140       if (constr_homology.gt.0) then
 141 c       write (iout,*) "About to call read_constr_homology"
 142 c       call flush(iout)
 143         call read_constr_homology
 144        write (iout,*) "Exit read_constr_homology"
 145        call flush(iout)
 146 cref        if (indpdb.gt.0 .or. pdbref) then
 147 cref          do i=1,2*nres
 148 cref            do j=1,3
 149 cref              c(j,i)=crefjlee(j,i)
 150 cref              cref(j,i)=crefjlee(j,i)
 151 cref            enddo
 152 cref          enddo
 153 cref        endif
 154       else
 155         homol_nset=0
 156       endif
 157
 158
 159       call setup_var
 160       call init_int_table
 161       if (ns.gt.0) then
 162         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 163      &  ' disulfide-bridging cysteines.'
 164         write (iout,'(20i4)') (iss(i),i=1,ns)
 165        if (dyn_ss) then
 166           write(iout,*)"Running with dynamic disulfide-bond formation"
 167        else
 168         write (iout,'(/a/)') 'Pre-formed links are:'
 169         do i=1,nss
 170           i1=ihpb(i)-nres
 171           i2=jhpb(i)-nres
 172           it1=itype(i1)
 173           it2=itype(i2)
 174          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 175      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 176      &    dhpb(i),ebr,forcon(i)
 177         enddo
 178       endif
 179       endif
 180       write (iout,'(a)')
 181       if (ns.gt.0.and.dyn_ss) then
 182           do i=nss+1,nhpb
 183             ihpb(i-nss)=ihpb(i)
 184             jhpb(i-nss)=jhpb(i)
 185             forcon(i-nss)=forcon(i)
 186             dhpb(i-nss)=dhpb(i)
 187           enddo
 188           nhpb=nhpb-nss
 189           nss=0
 190           call hpb_partition
 191           do i=1,ns
 192             dyn_ss_mask(iss(i))=.true.
 193           enddo
 194       endif
 195       return
 196       end
 197 c-----------------------------------------------------------------------------
 198       logical function seq_comp(itypea,itypeb,length)
 199       implicit none
 200       integer length,itypea(length),itypeb(length)
 201       integer i
 202       do i=1,length
 203         if (itypea(i).ne.itypeb(i)) then
 204           seq_comp=.false.
 205           return
 206         endif
 207       enddo
 208       seq_comp=.true.
 209       return
 210       end
 211 c-----------------------------------------------------------------------------
 212       subroutine read_bridge
 213 C Read information about disulfide bridges.
 214       implicit real*8 (a-h,o-z)
 215       include 'DIMENSIONS'
 216       include 'DIMENSIONS.ZSCOPT'
 217       include 'COMMON.IOUNITS'
 218       include 'COMMON.GEO'
 219       include 'COMMON.VAR'
 220       include 'COMMON.INTERACT'
 221       include 'COMMON.NAMES'
 222       include 'COMMON.CHAIN'
 223       include 'COMMON.FFIELD'
 224       include 'COMMON.SBRIDGE'
 225 C Read bridging residues.
 226       read (inp,*) ns,(iss(i),i=1,ns)
 227       print *,'ns=',ns
 228       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 229 C Check whether the specified bridging residues are cystines.
 230       do i=1,ns
 231         if (itype(iss(i)).ne.1) then
 232           write (iout,'(2a,i3,a)')
 233      &   'Do you REALLY think that the residue ',
 234      &    restyp(itype(iss(i))),i,
 235      &   ' can form a disulfide bridge?!!!'
 236           write (*,'(2a,i3,a)')
 237      &   'Do you REALLY think that the residue ',
 238      &    restyp(itype(iss(i))),i,
 239      &   ' can form a disulfide bridge?!!!'
 240          stop
 241         endif
 242       enddo
 243 C Read preformed bridges.
 244       if (ns.gt.0) then
 245       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 246       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 247       if (nss.gt.0) then
 248         nhpb=nss
 249 C Check if the residues involved in bridges are in the specified list of
 250 C bridging residues.
 251         do i=1,nss
 252           do j=1,i-1
 253             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 254      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 255               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 256      &      ' contains residues present in other pairs.'
 257               write (*,'(a,i3,a)') 'Disulfide pair',i,
 258      &      ' contains residues present in other pairs.'
 259               stop
 260             endif
 261           enddo
 262           do j=1,ns
 263             if (ihpb(i).eq.iss(j)) goto 10
 264           enddo
 265           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 266    10     continue
 267           do j=1,ns
 268             if (jhpb(i).eq.iss(j)) goto 20
 269           enddo
 270           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 271    20     continue
 272           dhpb(i)=dbr
 273           forcon(i)=fbr
 274         enddo
 275         do i=1,nss
 276           ihpb(i)=ihpb(i)+nres
 277           jhpb(i)=jhpb(i)+nres
 278         enddo
 279       endif
 280       endif
 281       return
 282       end
 283 c------------------------------------------------------------------------------
 284       subroutine read_angles(kanal,iscor,energ,iprot,*)
 285       implicit real*8 (a-h,o-z)
 286       include 'DIMENSIONS'
 287       include 'DIMENSIONS.ZSCOPT'
 288       include 'COMMON.INTERACT'
 289       include 'COMMON.SBRIDGE'
 290       include 'COMMON.GEO'
 291       include 'COMMON.VAR'
 292       include 'COMMON.CHAIN'
 293       include 'COMMON.IOUNITS'
 294       character*80 lineh
 295       read(kanal,'(a80)',end=10,err=10) lineh
 296       read(lineh(:5),*,err=8) ic
 297       read(lineh(6:),*,err=8) energ
 298       goto 9
 299     8 ic=1
 300       print *,'error, assuming e=1d10',lineh
 301       energ=1d10
 302       nss=0
 303     9 continue
 304       read(lineh(18:),*,end=10,err=10) nss
 305       IF (NSS.LT.9) THEN
 306         read (lineh(20:),*,end=10,err=10)
 307      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 308       ELSE
 309         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 310         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 311      &    I=9,NSS),iscor
 312       ENDIF
 313 c      print *,"energy",energ," iscor",iscor
 314       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 315       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 316       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 317       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 318       do i=1,nres
 319         theta(i)=deg2rad*theta(i)
 320         phi(i)=deg2rad*phi(i)
 321         alph(i)=deg2rad*alph(i)
 322         omeg(i)=deg2rad*omeg(i)
 323       enddo
 324       return
 325    10 return1
 326       end
 327 c====-------------------------------------------------------------------
 328       subroutine read_constr_homology
 329
 330       include 'DIMENSIONS'
 331       include 'DIMENSIONS.ZSCOPT'
 332       include 'DIMENSIONS.FREE'
 333 #ifdef MPI
 334       include 'mpif.h'
 335 #endif
 336       include 'COMMON.SETUP'
 337       include 'COMMON.CONTROL'
 338       include 'COMMON.CHAIN'
 339       include 'COMMON.IOUNITS'
 340       include 'COMMON.GEO'
 341       include 'COMMON.INTERACT'
 342       include 'COMMON.HOMRESTR'
 343 c
 344 c For new homol impl
 345 c
 346       include 'COMMON.VAR'
 347 c     include 'include_unres/COMMON.VAR'
 348 c
 349
 350 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 351 c    &                 dist_cut
 352 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 353 c    &    sigma_odl_temp(maxres,maxres,max_template)
 354       character*2 kic2
 355       character*24 model_ki_dist, model_ki_angle
 356       character*500 controlcard
 357       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 358       integer idomain(max_template,maxres)
 359       logical lprn /.true./
 360       integer ilen
 361       external ilen
 362       logical unres_pdb
 363 c
 364 c     FP - Nov. 2014 Temporary specifications for new vars
 365 c
 366       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 367       double precision, dimension (max_template,maxres) :: rescore
 368       double precision, dimension (max_template,maxres) :: rescore2
 369       character*24 tpl_k_rescore
 370 c -----------------------------------------------------------------
 371 c Reading multiple PDB ref structures and calculation of retraints
 372 c not using pre-computed ones stored in files model_ki_{dist,angle}
 373 c FP (Nov., 2014)
 374 c -----------------------------------------------------------------
 375 c
 376 c
 377 c Alternative: reading from input
 378       call card_concat(controlcard,.true.)
 379       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 380       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 381       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 382       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 383       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 384       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 385       call readi(controlcard,"HOMOL_SET",homol_nset,1)
 386       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 387       call readi(controlcard,"IHSET",ihset,1)
 388       if (homol_nset.gt.1)then
 389          call card_concat(controlcard,.true.)
 390          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 391          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 392           write(iout,*) "iset homology_weight "
 393 c         do i=1,homol_nset
 394 c          write(iout,*) i,waga_homology(i)
 395 c         enddo
 396          endif
 397          iset=mod(kolor,homol_nset)+1
 398       else
 399        iset=1
 400        waga_homology(1)=1.0
 401       endif
 402 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 403
 404 cd      write (iout,*) "nnt",nnt," nct",nct
 405 cd      call flush(iout)
 406
 407
 408       lim_odl=0
 409       lim_dih=0
 410 c
 411 c  New
 412 c
 413       lim_theta=0
 414       lim_xx=0
 415 c
 416 c  Reading HM global scores (prob not required)
 417 c
 418       do i = nnt,nct
 419         do k=1,constr_homology
 420          idomain(k,i)=0
 421         enddo
 422       enddo
 423 c     open (4,file="HMscore")
 424 c     do k=1,constr_homology
 425 c       read (4,*,end=521) hmscore_tmp
 426 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 427 c       write(*,*) "Model", k, ":", hmscore(k)
 428 c     enddo
 429 c521  continue
 430
 431       ii=0
 432       do i = nnt,nct-2
 433         do j=i+2,nct
 434         ii=ii+1
 435         ii_in_use(ii)=0
 436         enddo
 437       enddo
 438 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 439
 440       do k=1,constr_homology
 441
 442         read(inp,'(a)') pdbfile
 443 c  Next stament causes error upon compilation (?)
 444 c       if(me.eq.king.or. .not. out1file)
 445 c         write (iout,'(2a)') 'PDB data will be read from file ',
 446 c    &   pdbfile(:ilen(pdbfile))
 447          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 448      &  pdbfile(:ilen(pdbfile))
 449         open(ipdbin,file=pdbfile,status='old',err=33)
 450         goto 34
 451   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 452      &  pdbfile(:ilen(pdbfile))
 453         stop
 454   34    continue
 455 c        print *,'Begin reading pdb data'
 456 c
 457 c Files containing res sim or local scores (former containing sigmas)
 458 c
 459
 460         write(kic2,'(bz,i2.2)') k
 461
 462         tpl_k_rescore="template"//kic2//".sco"
 463
 464         unres_pdb=.false.
 465         call readpdb
 466 cref        do i=1,2*nres
 467 cref          do j=1,3
 468 cref            crefjlee(j,i)=c(j,i)
 469 cref          enddo
 470 cref        enddo
 471 #ifdef DEBUG
 472         do i=1,nres
 473           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 474      &      (crefjlee(j,i+nres),j=1,3)
 475         enddo
 476 #endif
 477         write (iout,*) "read_constr_homology: after reading pdb file"
 478         call flush(iout)
 479
 480 c
 481 c     Distance restraints
 482 c
 483 c          ... --> odl(k,ii)
 484 C Copy the coordinates from reference coordinates (?)
 485         do i=1,2*nres
 486           do j=1,3
 487 c            c(j,i)=cref(j,i)
 488 c           write (iout,*) "c(",j,i,") =",c(j,i)
 489           enddo
 490         enddo
 491 c
 492 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 493 c
 494 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 495           open (ientin,file=tpl_k_rescore,status='old')
 496           if (nnt.gt.1) rescore(k,1)=0.0d0
 497           do irec=nnt,maxdim ! loop for reading res sim
 498             if (read2sigma) then
 499              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 500      &                                idomain_tmp
 501              i_tmp=i_tmp+nnt-1
 502              idomain(k,i_tmp)=idomain_tmp
 503              rescore(k,i_tmp)=rescore_tmp
 504              rescore2(k,i_tmp)=rescore2_tmp
 505             else
 506              idomain(k,irec)=1
 507              read (ientin,*,end=1401) rescore_tmp
 508
 509 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 510              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 511 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 512             endif
 513           enddo
 514  1401   continue
 515         close (ientin)
 516         if (waga_dist.ne.0.0d0) then
 517           ii=0
 518           do i = nnt,nct-2
 519             do j=i+2,nct
 520
 521               x12=c(1,i)-c(1,j)
 522               y12=c(2,i)-c(2,j)
 523               z12=c(3,i)-c(3,j)
 524               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 525 c              write (iout,*) k,i,j,distal,dist2_cut
 526
 527             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 528      &            .and. distal.le.dist2_cut ) then
 529
 530               ii=ii+1
 531               ii_in_use(ii)=1
 532               l_homo(k,ii)=.true.
 533
 534 c             write (iout,*) "k",k
 535 c             write (iout,*) "i",i," j",j," constr_homology",
 536 c    &                       constr_homology
 537               ires_homo(ii)=i
 538               jres_homo(ii)=j
 539               odl(k,ii)=distal
 540               if (read2sigma) then
 541                 sigma_odl(k,ii)=0
 542                 do ik=i,j
 543                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 544                 enddo
 545                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 546                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 547      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 548               else
 549                 if (odl(k,ii).le.dist_cut) then
 550                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 551                 else
 552 #ifdef OLDSIGMA
 553                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 554      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 555 #else
 556                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 557      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 558 #endif
 559                 endif
 560               endif
 561               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 562             else
 563               ii=ii+1
 564               l_homo(k,ii)=.false.
 565             endif
 566             enddo
 567           enddo
 568         lim_odl=ii
 569         endif
 570 c
 571 c     Theta, dihedral and SC retraints
 572 c
 573         if (waga_angle.gt.0.0d0) then
 574 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 575 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 576 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 577 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 578 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 579 c    &                            sigma_dih(k,i+nnt-1)
 580 c         enddo
 581 c1402   continue
 582 c         close (ientin)
 583           do i = nnt+3,nct
 584             if (idomain(k,i).eq.0) then
 585                sigma_dih(k,i)=0.0
 586                cycle
 587             endif
 588             dih(k,i)=phiref(i) ! right?
 589 c           read (ientin,*) sigma_dih(k,i) ! original variant
 590 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 591 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 592 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 593 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 594 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 595
 596             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 597      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 598 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 599 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 600 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 601 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 602 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 603 c   Instead of res sim other local measure of b/b str reliability possible
 604             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 605 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 606           enddo
 607           lim_dih=nct-nnt-2
 608         endif
 609
 610         if (waga_theta.gt.0.0d0) then
 611 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 612 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 613 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 614 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 615 c    &                              sigma_theta(k,i+nnt-1)
 616 c         enddo
 617 c1403   continue
 618 c         close (ientin)
 619
 620           do i = nnt+2,nct ! right? without parallel.
 621 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 622 c         do i=ithet_start,ithet_end ! with FG parallel.
 623              if (idomain(k,i).eq.0) then
 624               sigma_theta(k,i)=0.0
 625               cycle
 626              endif
 627              thetatpl(k,i)=thetaref(i)
 628 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 629 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 630 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 631 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 632 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 633              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 634      &                        rescore(k,i-2))/3.0
 635 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 636              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 637
 638 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 639 c                             rescore(k,i-2) !  right expression ?
 640 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 641           enddo
 642         endif
 643         lim_theta=nct-nnt-1
 644
 645         if (waga_d.gt.0.0d0) then
 646 c       open (ientin,file=tpl_k_sigma_d,status='old')
 647 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 648 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 649 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 650 c    &                          sigma_d(k,i+nnt-1)
 651 c         enddo
 652 c1404   continue
 653
 654           do i = nnt,nct ! right? without parallel.
 655 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 656 c         do i=loc_start,loc_end ! with FG parallel.
 657                if (itype(i).eq.10) cycle
 658                if (idomain(k,i).eq.0 ) then
 659                   sigma_d(k,i)=0.0
 660                   cycle
 661                endif
 662                xxtpl(k,i)=xxref(i)
 663                yytpl(k,i)=yyref(i)
 664                zztpl(k,i)=zzref(i)
 665 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 666 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 667 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 668 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 669                sigma_d(k,i)=rescore(k,i) !  right expression ?
 670                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 671
 672 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 673 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 674 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 675                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
 676           enddo
 677           lim_xx=nct-nnt+1
 678         endif
 679       enddo
 680 c
 681 c remove distance restraints not used in any model from the list
 682 c shift data in all arrays
 683 c
 684       if (waga_dist.ne.0.0d0) then
 685         ii=0
 686         do i=nnt,nct-2
 687          do j=i+2,nct
 688           ii=ii+1
 689           if (ii_in_use(ii).eq.0) then
 690              do ki=ii,lim_odl-1
 691               ires_homo(ki)=ires_homo(ki+1)
 692               jres_homo(ki)=jres_homo(ki+1)
 693               ii_in_use(ki)=ii_in_use(ki+1)
 694               do k=1,constr_homology
 695                odl(k,ki)=odl(k,ki+1)
 696                sigma_odl(k,ki)=sigma_odl(k,ki+1)
 697                l_homo(k,ki)=l_homo(k,ki+1)
 698               enddo
 699              enddo
 700              ii=ii-1
 701              lim_odl=lim_odl-1
 702           endif
 703          enddo
 704         enddo
 705       endif
 706       if (constr_homology.gt.0) call homology_partition
 707       if (constr_homology.gt.0) call init_int_table
 708 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 709 cd     & "lim_xx=",lim_xx
 710 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 711 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 712 c
 713 c Print restraints
 714 c
 715       if (.not.lprn) return
 716 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 717       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 718        write (iout,*) "Distance restraints from templates"
 719        do ii=1,lim_odl
 720        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
 721      &  ii,ires_homo(ii),jres_homo(ii),
 722      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
 723      &  ki=1,constr_homology)
 724        enddo
 725        write (iout,*) "Dihedral angle restraints from templates"
 726        do i=nnt+3,lim_dih
 727         write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
 728      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 729        enddo
 730        write (iout,*) "Virtual-bond angle restraints from templates"
 731        do i=nnt+2,lim_theta
 732         write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
 733      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 734        enddo
 735        write (iout,*) "SC restraints from templates"
 736        do i=nnt,lim_xx
 737         write(iout,'(i5,100(4f8.2,4x))') i,
 738      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 739      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 740        enddo
 741       endif
 742 c -----------------------------------------------------------------
 743       return
 744       end
 745 c----------------------------------------------------------------------