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