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 'include_unres/COMMON.NAMES'
 343       include 'COMMON.HOMRESTR'
 344 c
 345 c For new homol impl
 346 c
 347       include 'COMMON.VAR'
 348 c     include 'include_unres/COMMON.VAR'
 349 c
 350
 351 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 352 c    &                 dist_cut
 353 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 354 c    &    sigma_odl_temp(maxres,maxres,max_template)
 355       character*2 kic2
 356       character*24 model_ki_dist, model_ki_angle
 357       character*500 controlcard
 358       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 359       integer idomain(max_template,maxres)
 360       logical lprn /.true./
 361       integer ilen
 362       external ilen
 363       logical unres_pdb,liiflag
 364 c
 365 c     FP - Nov. 2014 Temporary specifications for new vars
 366 c
 367       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
 368      &                 rescore3_tmp
 369       double precision, dimension (max_template,maxres) :: rescore
 370       double precision, dimension (max_template,maxres) :: rescore2
 371       double precision, dimension (max_template,maxres) :: rescore3
 372       character*24 tpl_k_rescore
 373 c -----------------------------------------------------------------
 374 c Reading multiple PDB ref structures and calculation of retraints
 375 c not using pre-computed ones stored in files model_ki_{dist,angle}
 376 c FP (Nov., 2014)
 377 c -----------------------------------------------------------------
 378 c
 379 c
 380 c Alternative: reading from input
 381       call card_concat(controlcard,.true.)
 382       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 383       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 384       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 385       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 386       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 387       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 388       call readi(controlcard,"HOMOL_SET",homol_nset,1)
 389       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 390       call readi(controlcard,"IHSET",ihset,1)
 391       if (homol_nset.gt.1)then
 392          call card_concat(controlcard,.true.)
 393          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 394          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 395           write(iout,*) "iset homology_weight "
 396 c         do i=1,homol_nset
 397 c          write(iout,*) i,waga_homology(i)
 398 c         enddo
 399          endif
 400          iset=mod(kolor,homol_nset)+1
 401       else
 402        iset=1
 403        waga_homology(1)=1.0
 404       endif
 405 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 406
 407 cd      write (iout,*) "nnt",nnt," nct",nct
 408 cd      call flush(iout)
 409
 410
 411       lim_odl=0
 412       lim_dih=0
 413 c
 414 c  New
 415 c
 416 c
 417 c  Reading HM global scores (prob not required)
 418 c
 419       do i = nnt,nct
 420         do k=1,constr_homology
 421          idomain(k,i)=0
 422         enddo
 423       enddo
 424 c     open (4,file="HMscore")
 425 c     do k=1,constr_homology
 426 c       read (4,*,end=521) hmscore_tmp
 427 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 428 c       write(*,*) "Model", k, ":", hmscore(k)
 429 c     enddo
 430 c521  continue
 431
 432       ii=0
 433       do i = nnt,nct-2
 434         do j=i+2,nct
 435         ii=ii+1
 436         ii_in_use(ii)=0
 437         enddo
 438       enddo
 439 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 440
 441       do k=1,constr_homology
 442
 443         read(inp,'(a)') pdbfile
 444 c  Next stament causes error upon compilation (?)
 445 c       if(me.eq.king.or. .not. out1file)
 446 c         write (iout,'(2a)') 'PDB data will be read from file ',
 447 c    &   pdbfile(:ilen(pdbfile))
 448          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 449      &  pdbfile(:ilen(pdbfile))
 450         open(ipdbin,file=pdbfile,status='old',err=33)
 451         goto 34
 452   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 453      &  pdbfile(:ilen(pdbfile))
 454         stop
 455   34    continue
 456 c        print *,'Begin reading pdb data'
 457 c
 458 c Files containing res sim or local scores (former containing sigmas)
 459 c
 460
 461         write(kic2,'(bz,i2.2)') k
 462
 463         tpl_k_rescore="template"//kic2//".sco"
 464
 465         unres_pdb=.false.
 466         call readpdb_template(k)
 467 cref        do i=1,2*nres
 468 cref          do j=1,3
 469 cref            crefjlee(j,i)=c(j,i)
 470 cref          enddo
 471 cref        enddo
 472 #ifdef DEBUG
 473         do i=1,nres
 474           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 475      &      (crefjlee(j,i+nres),j=1,3)
 476         enddo
 477 #endif
 478         write (iout,*) "read_constr_homology: after reading pdb file"
 479         call flush(iout)
 480
 481 c
 482 c     Distance restraints
 483 c
 484 c          ... --> odl(k,ii)
 485 C Copy the coordinates from reference coordinates (?)
 486         do i=1,2*nres
 487           do j=1,3
 488 c            c(j,i)=cref(j,i)
 489 c           write (iout,*) "c(",j,i,") =",c(j,i)
 490           enddo
 491         enddo
 492 c
 493 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 494 c
 495 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 496           open (ientin,file=tpl_k_rescore,status='old')
 497           if (nnt.gt.1) rescore(k,1)=0.0d0
 498           do irec=nnt,nct ! loop for reading res sim
 499             if (read2sigma) then
 500              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 501      &                                rescore3_tmp,idomain_tmp
 502              i_tmp=i_tmp+nnt-1
 503              idomain(k,i_tmp)=idomain_tmp
 504              rescore(k,i_tmp)=rescore_tmp
 505              rescore2(k,i_tmp)=rescore2_tmp
 506              rescore3(k,i_tmp)=rescore3_tmp
 507              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
 508      &                      i_tmp,rescore2_tmp,rescore_tmp,
 509      &                                rescore3_tmp,idomain_tmp
 510             else
 511              idomain(k,irec)=1
 512              read (ientin,*,end=1401) rescore_tmp
 513
 514 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 515              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 516 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 517             endif
 518           enddo
 519  1401   continue
 520         close (ientin)
 521         if (waga_dist.ne.0.0d0) then
 522           ii=0
 523           do i = nnt,nct-2
 524             do j=i+2,nct
 525
 526               x12=c(1,i)-c(1,j)
 527               y12=c(2,i)-c(2,j)
 528               z12=c(3,i)-c(3,j)
 529               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 530 c              write (iout,*) k,i,j,distal,dist2_cut
 531
 532             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 533      &            .and. distal.le.dist2_cut ) then
 534
 535               ii=ii+1
 536               ii_in_use(ii)=1
 537               l_homo(k,ii)=.true.
 538
 539 c             write (iout,*) "k",k
 540 c             write (iout,*) "i",i," j",j," constr_homology",
 541 c    &                       constr_homology
 542               ires_homo(ii)=i
 543               jres_homo(ii)=j
 544               odl(k,ii)=distal
 545               if (read2sigma) then
 546                 sigma_odl(k,ii)=0
 547                 do ik=i,j
 548                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 549                 enddo
 550                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 551                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 552      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 553               else
 554                 if (odl(k,ii).le.dist_cut) then
 555                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 556                 else
 557 #ifdef OLDSIGMA
 558                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 559      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 560 #else
 561                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 562      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 563 #endif
 564                 endif
 565               endif
 566               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 567             else
 568               ii=ii+1
 569               l_homo(k,ii)=.false.
 570             endif
 571             enddo
 572           enddo
 573         lim_odl=ii
 574         endif
 575 c
 576 c     Theta, dihedral and SC retraints
 577 c
 578         if (waga_angle.gt.0.0d0) then
 579 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 580 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 581 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 582 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 583 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 584 c    &                            sigma_dih(k,i+nnt-1)
 585 c         enddo
 586 c1402   continue
 587 c         close (ientin)
 588           do i = nnt+3,nct
 589             if (idomain(k,i).eq.0) then
 590                sigma_dih(k,i)=0.0
 591                cycle
 592             endif
 593             dih(k,i)=phiref(i) ! right?
 594 c           read (ientin,*) sigma_dih(k,i) ! original variant
 595 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 596 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 597 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 598 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 599 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 600
 601             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 602      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 603 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 604 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 605 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 606 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 607 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 608 c   Instead of res sim other local measure of b/b str reliability possible
 609             if (sigma_dih(k,i).ne.0)
 610      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 611 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 612           enddo
 613           lim_dih=nct-nnt-2
 614         endif
 615
 616         if (waga_theta.gt.0.0d0) then
 617 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 618 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 619 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 620 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 621 c    &                              sigma_theta(k,i+nnt-1)
 622 c         enddo
 623 c1403   continue
 624 c         close (ientin)
 625
 626           do i = nnt+2,nct ! right? without parallel.
 627 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 628 c         do i=ithet_start,ithet_end ! with FG parallel.
 629              if (idomain(k,i).eq.0) then
 630               sigma_theta(k,i)=0.0
 631               cycle
 632              endif
 633              thetatpl(k,i)=thetaref(i)
 634 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 635 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 636 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 637 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 638 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 639              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 640      &                        rescore(k,i-2))/3.0
 641 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 642              if (sigma_theta(k,i).ne.0)
 643      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 644
 645 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 646 c                             rescore(k,i-2) !  right expression ?
 647 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 648           enddo
 649         endif
 650
 651         if (waga_d.gt.0.0d0) then
 652 c       open (ientin,file=tpl_k_sigma_d,status='old')
 653 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 654 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 655 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 656 c    &                          sigma_d(k,i+nnt-1)
 657 c         enddo
 658 c1404   continue
 659
 660           do i = nnt,nct ! right? without parallel.
 661 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 662 c         do i=loc_start,loc_end ! with FG parallel.
 663                if (itype(i).eq.10) cycle
 664                if (idomain(k,i).eq.0 ) then
 665                   sigma_d(k,i)=0.0
 666                   cycle
 667                endif
 668                xxtpl(k,i)=xxref(i)
 669                yytpl(k,i)=yyref(i)
 670                zztpl(k,i)=zzref(i)
 671 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 672 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 673 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 674 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 675                sigma_d(k,i)=rescore3(k,i) !  right expression ?
 676                if (sigma_d(k,i).ne.0)
 677      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 678
 679 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 680 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 681 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 682           enddo
 683         endif
 684       enddo
 685 c
 686 c remove distance restraints not used in any model from the list
 687 c shift data in all arrays
 688 c
 689       if (waga_dist.ne.0.0d0) then
 690         ii=0
 691         liiflag=.true.
 692         do i=nnt,nct-2
 693          do j=i+2,nct
 694           ii=ii+1
 695           if (ii_in_use(ii).eq.0.and.liiflag) then
 696             liiflag=.false.
 697             iistart=ii
 698           endif
 699           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
 700      &                   .not.liiflag.and.ii.eq.lim_odl) then
 701              if (ii.eq.lim_odl) then
 702               iishift=ii-iistart+1
 703              else
 704               iishift=ii-iistart
 705              endif
 706              liiflag=.true.
 707              do ki=iistart,lim_odl-iishift
 708               ires_homo(ki)=ires_homo(ki+iishift)
 709               jres_homo(ki)=jres_homo(ki+iishift)
 710               ii_in_use(ki)=ii_in_use(ki+iishift)
 711               do k=1,constr_homology
 712                odl(k,ki)=odl(k,ki+iishift)
 713                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
 714                l_homo(k,ki)=l_homo(k,ki+iishift)
 715               enddo
 716              enddo
 717              ii=ii-iishift
 718              lim_odl=lim_odl-iishift
 719           endif
 720          enddo
 721         enddo
 722       endif
 723       if (constr_homology.gt.0) call homology_partition
 724       if (constr_homology.gt.0) call init_int_table
 725 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 726 cd     & "lim_xx=",lim_xx
 727 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 728 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 729 c
 730 c Print restraints
 731 c
 732       if (.not.lprn) return
 733 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 734       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 735        write (iout,*) "Distance restraints from templates"
 736        do ii=1,lim_odl
 737        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
 738      &  ii,ires_homo(ii),jres_homo(ii),
 739      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
 740      &  ki=1,constr_homology)
 741        enddo
 742        write (iout,*) "Dihedral angle restraints from templates"
 743        do i=nnt+3,nct
 744         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 745      &      (rad2deg*dih(ki,i),
 746      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 747        enddo
 748        write (iout,*) "Virtual-bond angle restraints from templates"
 749        do i=nnt+2,nct
 750         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 751      &      (rad2deg*thetatpl(ki,i),
 752      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 753        enddo
 754        write (iout,*) "SC restraints from templates"
 755        do i=nnt,nct
 756         write(iout,'(i5,100(4f8.2,4x))') i,
 757      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 758      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 759        enddo
 760       endif
 761 c -----------------------------------------------------------------
 762       return
 763       end
 764 c----------------------------------------------------------------------