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       write (iout,*) "calling read_saxs_consrtr",nsaxs
 141       if (nsaxs.gt.0) call read_saxs_constr
 142
 143       if (constr_homology.gt.0) then
 144 c       write (iout,*) "About to call read_constr_homology"
 145 c       call flush(iout)
 146         call read_constr_homology
 147        write (iout,*) "Exit read_constr_homology"
 148        call flush(iout)
 149 cref        if (indpdb.gt.0 .or. pdbref) then
 150 cref          do i=1,2*nres
 151 cref            do j=1,3
 152 cref              c(j,i)=crefjlee(j,i)
 153 cref              cref(j,i)=crefjlee(j,i)
 154 cref            enddo
 155 cref          enddo
 156 cref        endif
 157       else
 158         homol_nset=0
 159       endif
 160
 161
 162       call setup_var
 163       call init_int_table
 164       if (ns.gt.0) then
 165         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 166      &  ' disulfide-bridging cysteines.'
 167         write (iout,'(20i4)') (iss(i),i=1,ns)
 168        if (dyn_ss) then
 169           write(iout,*)"Running with dynamic disulfide-bond formation"
 170        else
 171         write (iout,'(/a/)') 'Pre-formed links are:'
 172         do i=1,nss
 173           i1=ihpb(i)-nres
 174           i2=jhpb(i)-nres
 175           it1=itype(i1)
 176           it2=itype(i2)
 177          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 178      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 179      &    dhpb(i),ebr,forcon(i)
 180         enddo
 181       endif
 182       endif
 183       write (iout,'(a)')
 184       if (ns.gt.0.and.dyn_ss) then
 185           do i=nss+1,nhpb
 186             ihpb(i-nss)=ihpb(i)
 187             jhpb(i-nss)=jhpb(i)
 188             forcon(i-nss)=forcon(i)
 189             dhpb(i-nss)=dhpb(i)
 190           enddo
 191           nhpb=nhpb-nss
 192           nss=0
 193           call hpb_partition
 194           do i=1,ns
 195             dyn_ss_mask(iss(i))=.true.
 196           enddo
 197       endif
 198       return
 199       end
 200 c-----------------------------------------------------------------------------
 201       logical function seq_comp(itypea,itypeb,length)
 202       implicit none
 203       integer length,itypea(length),itypeb(length)
 204       integer i
 205       do i=1,length
 206         if (itypea(i).ne.itypeb(i)) then
 207           seq_comp=.false.
 208           return
 209         endif
 210       enddo
 211       seq_comp=.true.
 212       return
 213       end
 214 c-----------------------------------------------------------------------------
 215       subroutine read_bridge
 216 C Read information about disulfide bridges.
 217       implicit real*8 (a-h,o-z)
 218       include 'DIMENSIONS'
 219       include 'DIMENSIONS.ZSCOPT'
 220       include 'COMMON.IOUNITS'
 221       include 'COMMON.GEO'
 222       include 'COMMON.VAR'
 223       include 'COMMON.INTERACT'
 224       include 'COMMON.NAMES'
 225       include 'COMMON.CHAIN'
 226       include 'COMMON.FFIELD'
 227       include 'COMMON.SBRIDGE'
 228 C Read bridging residues.
 229       read (inp,*) ns,(iss(i),i=1,ns)
 230       print *,'ns=',ns
 231       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 232 C Check whether the specified bridging residues are cystines.
 233       do i=1,ns
 234         if (itype(iss(i)).ne.1) then
 235           write (iout,'(2a,i3,a)')
 236      &   'Do you REALLY think that the residue ',
 237      &    restyp(itype(iss(i))),i,
 238      &   ' can form a disulfide bridge?!!!'
 239           write (*,'(2a,i3,a)')
 240      &   'Do you REALLY think that the residue ',
 241      &    restyp(itype(iss(i))),i,
 242      &   ' can form a disulfide bridge?!!!'
 243          stop
 244         endif
 245       enddo
 246 C Read preformed bridges.
 247       if (ns.gt.0) then
 248       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 249       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 250       if (nss.gt.0) then
 251         nhpb=nss
 252 C Check if the residues involved in bridges are in the specified list of
 253 C bridging residues.
 254         do i=1,nss
 255           do j=1,i-1
 256             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 257      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 258               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 259      &      ' contains residues present in other pairs.'
 260               write (*,'(a,i3,a)') 'Disulfide pair',i,
 261      &      ' contains residues present in other pairs.'
 262               stop
 263             endif
 264           enddo
 265           do j=1,ns
 266             if (ihpb(i).eq.iss(j)) goto 10
 267           enddo
 268           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 269    10     continue
 270           do j=1,ns
 271             if (jhpb(i).eq.iss(j)) goto 20
 272           enddo
 273           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 274    20     continue
 275           dhpb(i)=dbr
 276           forcon(i)=fbr
 277         enddo
 278         do i=1,nss
 279           ihpb(i)=ihpb(i)+nres
 280           jhpb(i)=jhpb(i)+nres
 281         enddo
 282       endif
 283       endif
 284       return
 285       end
 286 c------------------------------------------------------------------------------
 287       subroutine read_angles(kanal,iscor,energ,iprot,*)
 288       implicit real*8 (a-h,o-z)
 289       include 'DIMENSIONS'
 290       include 'DIMENSIONS.ZSCOPT'
 291       include 'COMMON.INTERACT'
 292       include 'COMMON.SBRIDGE'
 293       include 'COMMON.GEO'
 294       include 'COMMON.VAR'
 295       include 'COMMON.CHAIN'
 296       include 'COMMON.IOUNITS'
 297       character*80 lineh
 298       read(kanal,'(a80)',end=10,err=10) lineh
 299       read(lineh(:5),*,err=8) ic
 300       read(lineh(6:),*,err=8) energ
 301       goto 9
 302     8 ic=1
 303       print *,'error, assuming e=1d10',lineh
 304       energ=1d10
 305       nss=0
 306     9 continue
 307       read(lineh(18:),*,end=10,err=10) nss
 308       IF (NSS.LT.9) THEN
 309         read (lineh(20:),*,end=10,err=10)
 310      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 311       ELSE
 312         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 313         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 314      &    I=9,NSS),iscor
 315       ENDIF
 316 c      print *,"energy",energ," iscor",iscor
 317       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 318       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 319       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 320       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 321       do i=1,nres
 322         theta(i)=deg2rad*theta(i)
 323         phi(i)=deg2rad*phi(i)
 324         alph(i)=deg2rad*alph(i)
 325         omeg(i)=deg2rad*omeg(i)
 326       enddo
 327       return
 328    10 return1
 329       end
 330 c-------------------------------------------------------------------------------
 331       subroutine read_saxs_constr
 332       implicit real*8 (a-h,o-z)
 333       include 'DIMENSIONS'
 334       include 'DIMENSIONS.ZSCOPT'
 335       include 'DIMENSIONS.FREE'
 336 #ifdef MPI
 337       include 'mpif.h'
 338 #endif
 339       include 'COMMON.CONTROL'
 340       include 'COMMON.CHAIN'
 341       include 'COMMON.IOUNITS'
 342       include 'COMMON.SBRIDGE'
 343       double precision cm(3)
 344 c      read(inp,*) nsaxs
 345       write (iout,*) "Calling read_saxs nsaxs",nsaxs
 346       call flush(iout)
 347       if (saxs_mode.eq.0) then
 348 c SAXS distance distribution
 349       do i=1,nsaxs
 350         read(inp,*) distsaxs(i),Psaxs(i)
 351       enddo
 352       Cnorm = 0.0d0
 353       do i=1,nsaxs
 354         Cnorm = Cnorm + Psaxs(i)
 355       enddo
 356       write (iout,*) "Cnorm",Cnorm
 357       do i=1,nsaxs
 358         Psaxs(i)=Psaxs(i)/Cnorm
 359       enddo
 360       write (iout,*) "Normalized distance distribution from SAXS"
 361       do i=1,nsaxs
 362         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
 363       enddo
 364       else
 365 c SAXS "spheres".
 366       do i=1,nsaxs
 367         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
 368       enddo
 369       do j=1,3
 370         cm(j)=0.0d0
 371       enddo
 372       do i=1,nsaxs
 373         do j=1,3
 374           cm(j)=cm(j)+Csaxs(j,i)
 375         enddo
 376       enddo
 377       do j=1,3
 378         cm(j)=cm(j)/nsaxs
 379       enddo
 380       do i=1,nsaxs
 381         do j=1,3
 382           Csaxs(j,i)=Csaxs(j,i)-cm(j)
 383         enddo
 384       enddo
 385       write (iout,*) "SAXS sphere coordinates"
 386       do i=1,nsaxs
 387         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
 388       enddo
 389       endif
 390       return
 391       end
 392 c====-------------------------------------------------------------------
 393       subroutine read_constr_homology
 394
 395       include 'DIMENSIONS'
 396       include 'DIMENSIONS.ZSCOPT'
 397       include 'DIMENSIONS.FREE'
 398 #ifdef MPI
 399       include 'mpif.h'
 400 #endif
 401       include 'COMMON.SETUP'
 402       include 'COMMON.CONTROL'
 403       include 'COMMON.CHAIN'
 404       include 'COMMON.IOUNITS'
 405       include 'COMMON.GEO'
 406       include 'COMMON.INTERACT'
 407       include 'include_unres/COMMON.NAMES'
 408       include 'COMMON.HOMRESTR'
 409 c
 410 c For new homol impl
 411 c
 412       include 'COMMON.VAR'
 413 c     include 'include_unres/COMMON.VAR'
 414 c
 415
 416 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 417 c    &                 dist_cut
 418 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 419 c    &    sigma_odl_temp(maxres,maxres,max_template)
 420       character*2 kic2
 421       character*24 model_ki_dist, model_ki_angle
 422       character*500 controlcard
 423       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 424       integer idomain(max_template,maxres)
 425       logical lprn /.true./
 426       integer ilen
 427       external ilen
 428       logical unres_pdb,liiflag
 429 c
 430 c     FP - Nov. 2014 Temporary specifications for new vars
 431 c
 432       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
 433      &                 rescore3_tmp
 434       double precision, dimension (max_template,maxres) :: rescore
 435       double precision, dimension (max_template,maxres) :: rescore2
 436       double precision, dimension (max_template,maxres) :: rescore3
 437       character*24 tpl_k_rescore
 438 c -----------------------------------------------------------------
 439 c Reading multiple PDB ref structures and calculation of retraints
 440 c not using pre-computed ones stored in files model_ki_{dist,angle}
 441 c FP (Nov., 2014)
 442 c -----------------------------------------------------------------
 443 c
 444 c
 445 c Alternative: reading from input
 446       call card_concat(controlcard,.true.)
 447       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 448       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 449       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 450       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 451       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 452       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 453       call readi(controlcard,"HOMOL_SET",homol_nset,1)
 454       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 455       call readi(controlcard,"IHSET",ihset,1)
 456       if (homol_nset.gt.1)then
 457          call card_concat(controlcard,.true.)
 458          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 459          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 460           write(iout,*) "iset homology_weight "
 461 c         do i=1,homol_nset
 462 c          write(iout,*) i,waga_homology(i)
 463 c         enddo
 464          endif
 465          iset=mod(kolor,homol_nset)+1
 466       else
 467        iset=1
 468        waga_homology(1)=1.0
 469       endif
 470 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 471
 472 cd      write (iout,*) "nnt",nnt," nct",nct
 473 cd      call flush(iout)
 474
 475
 476       lim_odl=0
 477       lim_dih=0
 478 c
 479 c  New
 480 c
 481 c
 482 c  Reading HM global scores (prob not required)
 483 c
 484       do i = nnt,nct
 485         do k=1,constr_homology
 486          idomain(k,i)=0
 487         enddo
 488       enddo
 489 c     open (4,file="HMscore")
 490 c     do k=1,constr_homology
 491 c       read (4,*,end=521) hmscore_tmp
 492 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 493 c       write(*,*) "Model", k, ":", hmscore(k)
 494 c     enddo
 495 c521  continue
 496
 497       ii=0
 498       do i = nnt,nct-2
 499         do j=i+2,nct
 500         ii=ii+1
 501         ii_in_use(ii)=0
 502         enddo
 503       enddo
 504 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 505
 506       do k=1,constr_homology
 507
 508         read(inp,'(a)') pdbfile
 509 c  Next stament causes error upon compilation (?)
 510 c       if(me.eq.king.or. .not. out1file)
 511 c         write (iout,'(2a)') 'PDB data will be read from file ',
 512 c    &   pdbfile(:ilen(pdbfile))
 513          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 514      &  pdbfile(:ilen(pdbfile))
 515         open(ipdbin,file=pdbfile,status='old',err=33)
 516         goto 34
 517   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 518      &  pdbfile(:ilen(pdbfile))
 519         stop
 520   34    continue
 521 c        print *,'Begin reading pdb data'
 522 c
 523 c Files containing res sim or local scores (former containing sigmas)
 524 c
 525
 526         write(kic2,'(bz,i2.2)') k
 527
 528         tpl_k_rescore="template"//kic2//".sco"
 529
 530         unres_pdb=.false.
 531         call readpdb_template(k)
 532 cref        do i=1,2*nres
 533 cref          do j=1,3
 534 cref            crefjlee(j,i)=c(j,i)
 535 cref          enddo
 536 cref        enddo
 537 #ifdef DEBUG
 538         do i=1,nres
 539           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 540      &      (crefjlee(j,i+nres),j=1,3)
 541         enddo
 542 #endif
 543         write (iout,*) "read_constr_homology: after reading pdb file"
 544         call flush(iout)
 545
 546 c
 547 c     Distance restraints
 548 c
 549 c          ... --> odl(k,ii)
 550 C Copy the coordinates from reference coordinates (?)
 551         do i=1,2*nres
 552           do j=1,3
 553 c            c(j,i)=cref(j,i)
 554 c           write (iout,*) "c(",j,i,") =",c(j,i)
 555           enddo
 556         enddo
 557 c
 558 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 559 c
 560 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 561           open (ientin,file=tpl_k_rescore,status='old')
 562           if (nnt.gt.1) rescore(k,1)=0.0d0
 563           do irec=nnt,nct ! loop for reading res sim
 564             if (read2sigma) then
 565              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 566      &                                rescore3_tmp,idomain_tmp
 567              i_tmp=i_tmp+nnt-1
 568              idomain(k,i_tmp)=idomain_tmp
 569              rescore(k,i_tmp)=rescore_tmp
 570              rescore2(k,i_tmp)=rescore2_tmp
 571              rescore3(k,i_tmp)=rescore3_tmp
 572              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
 573      &                      i_tmp,rescore2_tmp,rescore_tmp,
 574      &                                rescore3_tmp,idomain_tmp
 575             else
 576              idomain(k,irec)=1
 577              read (ientin,*,end=1401) rescore_tmp
 578
 579 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 580              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 581 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 582             endif
 583           enddo
 584  1401   continue
 585         close (ientin)
 586         if (waga_dist.ne.0.0d0) then
 587           ii=0
 588           do i = nnt,nct-2
 589             do j=i+2,nct
 590
 591               x12=c(1,i)-c(1,j)
 592               y12=c(2,i)-c(2,j)
 593               z12=c(3,i)-c(3,j)
 594               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 595 c              write (iout,*) k,i,j,distal,dist2_cut
 596
 597             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 598      &            .and. distal.le.dist2_cut ) then
 599
 600               ii=ii+1
 601               ii_in_use(ii)=1
 602               l_homo(k,ii)=.true.
 603
 604 c             write (iout,*) "k",k
 605 c             write (iout,*) "i",i," j",j," constr_homology",
 606 c    &                       constr_homology
 607               ires_homo(ii)=i
 608               jres_homo(ii)=j
 609               odl(k,ii)=distal
 610               if (read2sigma) then
 611                 sigma_odl(k,ii)=0
 612                 do ik=i,j
 613                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 614                 enddo
 615                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 616                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 617      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 618               else
 619                 if (odl(k,ii).le.dist_cut) then
 620                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 621                 else
 622 #ifdef OLDSIGMA
 623                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 624      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 625 #else
 626                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 627      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 628 #endif
 629                 endif
 630               endif
 631               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 632             else
 633               ii=ii+1
 634               l_homo(k,ii)=.false.
 635             endif
 636             enddo
 637           enddo
 638         lim_odl=ii
 639         endif
 640 c
 641 c     Theta, dihedral and SC retraints
 642 c
 643         if (waga_angle.gt.0.0d0) then
 644 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 645 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 646 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 647 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 648 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 649 c    &                            sigma_dih(k,i+nnt-1)
 650 c         enddo
 651 c1402   continue
 652 c         close (ientin)
 653           do i = nnt+3,nct
 654             if (idomain(k,i).eq.0) then
 655                sigma_dih(k,i)=0.0
 656                cycle
 657             endif
 658             dih(k,i)=phiref(i) ! right?
 659 c           read (ientin,*) sigma_dih(k,i) ! original variant
 660 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 661 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 662 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 663 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 664 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 665
 666             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 667      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 668 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 669 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 670 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 671 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 672 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 673 c   Instead of res sim other local measure of b/b str reliability possible
 674             if (sigma_dih(k,i).ne.0)
 675      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 676 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 677           enddo
 678           lim_dih=nct-nnt-2
 679         endif
 680
 681         if (waga_theta.gt.0.0d0) then
 682 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 683 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 684 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 685 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 686 c    &                              sigma_theta(k,i+nnt-1)
 687 c         enddo
 688 c1403   continue
 689 c         close (ientin)
 690
 691           do i = nnt+2,nct ! right? without parallel.
 692 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 693 c         do i=ithet_start,ithet_end ! with FG parallel.
 694              if (idomain(k,i).eq.0) then
 695               sigma_theta(k,i)=0.0
 696               cycle
 697              endif
 698              thetatpl(k,i)=thetaref(i)
 699 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 700 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 701 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 702 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 703 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 704              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 705      &                        rescore(k,i-2))/3.0
 706 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 707              if (sigma_theta(k,i).ne.0)
 708      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 709
 710 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 711 c                             rescore(k,i-2) !  right expression ?
 712 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 713           enddo
 714         endif
 715
 716         if (waga_d.gt.0.0d0) then
 717 c       open (ientin,file=tpl_k_sigma_d,status='old')
 718 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 719 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 720 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 721 c    &                          sigma_d(k,i+nnt-1)
 722 c         enddo
 723 c1404   continue
 724
 725           do i = nnt,nct ! right? without parallel.
 726 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 727 c         do i=loc_start,loc_end ! with FG parallel.
 728                if (itype(i).eq.10) cycle
 729                if (idomain(k,i).eq.0 ) then
 730                   sigma_d(k,i)=0.0
 731                   cycle
 732                endif
 733                xxtpl(k,i)=xxref(i)
 734                yytpl(k,i)=yyref(i)
 735                zztpl(k,i)=zzref(i)
 736 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 737 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 738 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 739 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 740                sigma_d(k,i)=rescore3(k,i) !  right expression ?
 741                if (sigma_d(k,i).ne.0)
 742      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 743
 744 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 745 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 746 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 747           enddo
 748         endif
 749       enddo
 750 c
 751 c remove distance restraints not used in any model from the list
 752 c shift data in all arrays
 753 c
 754       if (waga_dist.ne.0.0d0) then
 755         ii=0
 756         liiflag=.true.
 757         do i=nnt,nct-2
 758          do j=i+2,nct
 759           ii=ii+1
 760           if (ii_in_use(ii).eq.0.and.liiflag) then
 761             liiflag=.false.
 762             iistart=ii
 763           endif
 764           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
 765      &                   .not.liiflag.and.ii.eq.lim_odl) then
 766              if (ii.eq.lim_odl) then
 767               iishift=ii-iistart+1
 768              else
 769               iishift=ii-iistart
 770              endif
 771              liiflag=.true.
 772              do ki=iistart,lim_odl-iishift
 773               ires_homo(ki)=ires_homo(ki+iishift)
 774               jres_homo(ki)=jres_homo(ki+iishift)
 775               ii_in_use(ki)=ii_in_use(ki+iishift)
 776               do k=1,constr_homology
 777                odl(k,ki)=odl(k,ki+iishift)
 778                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
 779                l_homo(k,ki)=l_homo(k,ki+iishift)
 780               enddo
 781              enddo
 782              ii=ii-iishift
 783              lim_odl=lim_odl-iishift
 784           endif
 785          enddo
 786         enddo
 787       endif
 788       if (constr_homology.gt.0) call homology_partition
 789       if (constr_homology.gt.0) call init_int_table
 790 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 791 cd     & "lim_xx=",lim_xx
 792 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 793 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 794 c
 795 c Print restraints
 796 c
 797       if (.not.lprn) return
 798 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 799       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 800        write (iout,*) "Distance restraints from templates"
 801        do ii=1,lim_odl
 802        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
 803      &  ii,ires_homo(ii),jres_homo(ii),
 804      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
 805      &  ki=1,constr_homology)
 806        enddo
 807        write (iout,*) "Dihedral angle restraints from templates"
 808        do i=nnt+3,nct
 809         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 810      &      (rad2deg*dih(ki,i),
 811      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 812        enddo
 813        write (iout,*) "Virtual-bond angle restraints from templates"
 814        do i=nnt+2,nct
 815         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 816      &      (rad2deg*thetatpl(ki,i),
 817      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 818        enddo
 819        write (iout,*) "SC restraints from templates"
 820        do i=nnt,nct
 821         write(iout,'(i5,100(4f8.2,4x))') i,
 822      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 823      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 824        enddo
 825       endif
 826 c -----------------------------------------------------------------
 827       return
 828       end
 829 c----------------------------------------------------------------------