source/wham/src/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 'COMMON.IOUNITS'
   9       include 'COMMON.GEO'
  10       include 'COMMON.VAR'
  11 c     include 'include_unres/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 C     Bartek
  32       call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
  33       call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
  34       call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
  35       call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
  36       write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
  37      &  " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
  38       r0_corr=cutoff_corr-delt_corr
  39       call readi(controlcard,"NRES",nres,0)
  40       iscode=index(controlcard,"ONE_LETTER")
  41       if (nres.le.0) then
  42         write (iout,*) "Error: no residues in molecule"
  43         return1
  44       endif
  45       if (nres.gt.maxres) then
  46         write (iout,*) "Error: too many residues",nres,maxres
  47       endif
  48       write(iout,*) 'nres=',nres
  49 C Read sequence of the protein
  50       if (iscode.gt.0) then
  51         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
  52       else
  53         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
  54       endif
  55 C Convert sequence to numeric code
  56       do i=1,nres
  57         itype(i)=rescode(i,sequence(i),iscode)
  58       enddo
  59       write (iout,*) "Numeric code:"
  60       write (iout,'(20i4)') (itype(i),i=1,nres)
  61       do i=1,nres-1
  62 #ifdef PROCOR
  63         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
  64 #else
  65         if (itype(i).eq.21) then
  66 #endif
  67           itel(i)=0
  68 #ifdef PROCOR
  69         else if (itype(i+1).ne.20) then
  70 #else
  71         else if (itype(i).ne.20) then
  72 #endif
  73           itel(i)=1
  74         else
  75           itel(i)=2
  76         endif
  77       enddo
  78       call read_bridge
  79
  80       if (with_dihed_constr) then
  81
  82       read (inp,*) ndih_constr
  83       if (ndih_constr.gt.0) then
  84         read (inp,*) ftors
  85         write (iout,*) 'FTORS',ftors
  86         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
  87         write (iout,*)
  88      &   'There are',ndih_constr,' constraints on phi angles.'
  89         do i=1,ndih_constr
  90           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
  91         enddo
  92         do i=1,ndih_constr
  93           phi0(i)=deg2rad*phi0(i)
  94           drange(i)=deg2rad*drange(i)
  95         enddo
  96       endif
  97
  98       endif
  99
 100       nnt=1
 101       nct=nres
 102       if (itype(1).eq.21) nnt=2
 103       if (itype(nres).eq.21) nct=nct-1
 104       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 105
 106 C     Juyong:READ init_vars
 107 C     Initialize variables!
 108 C     Juyong:READ read_info
 109 C     READ fragment information!!
 110 C     both routines should be in dfa.F file!!
 111
 112       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 113      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 114        write (iout,*) "Calling init_dfa_vars"
 115        call flush(iout)
 116        call init_dfa_vars
 117        write (iout,*) 'init_dfa_vars finished!'
 118        call flush(iout)
 119        call read_dfa_info
 120        write (iout,*) 'read_dfa_info finished!'
 121        call flush(iout)
 122       endif
 123
 124 c Read distance restraints
 125       if (constr_dist.gt.0) then
 126         if (refstr) call read_ref_structure(*11)
 127         call read_dist_constr
 128         call hpb_partition
 129       endif
 130
 131       if (constr_homology.gt.0) then
 132 c       write (iout,*) "About to call read_constr_homology"
 133 c       call flush(iout)
 134         call read_constr_homology
 135 c       write (iout,*) "Exit read_constr_homology"
 136 c       call flush(iout)
 137         if (indpdb.gt.0 .or. pdbref) then
 138           do i=1,2*nres
 139             do j=1,3
 140               c(j,i)=crefjlee(j,i)
 141               cref(j,i)=crefjlee(j,i)
 142             enddo
 143           enddo
 144         endif
 145 #ifdef DEBUG
 146         write (iout,*) "Array C"
 147         do i=1,nres
 148           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
 149      &      (c(j,i+nres),j=1,3)
 150         enddo
 151         write (iout,*) "Array Cref"
 152         do i=1,nres
 153           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
 154      &      (cref(j,i+nres),j=1,3)
 155         enddo
 156 #endif
 157 #ifdef DEBUG
 158        call int_from_cart1(.false.)
 159        call sc_loc_geom(.false.)
 160        do i=1,nres
 161          thetaref(i)=theta(i)
 162          phiref(i)=phi(i)
 163          write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
 164        enddo
 165        do i=1,nres-1
 166          do j=1,3
 167            dc(j,i)=c(j,i+1)-c(j,i)
 168            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 169          enddo
 170        enddo
 171        do i=2,nres-1
 172          do j=1,3
 173            dc(j,i+nres)=c(j,i+nres)-c(j,i)
 174            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 175          enddo
 176        enddo
 177 #endif
 178       else
 179         homol_nset=0
 180       endif
 181
 182
 183       call setup_var
 184       call init_int_table
 185       if (ns.gt.0) then
 186         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 187      &  ' disulfide-bridging cysteines.'
 188         write (iout,'(20i4)') (iss(i),i=1,ns)
 189         write (iout,'(/a/)') 'Pre-formed links are:'
 190         do i=1,nss
 191           i1=ihpb(i)-nres
 192           i2=jhpb(i)-nres
 193           it1=itype(i1)
 194           it2=itype(i2)
 195          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 196      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 197      &    dhpb(i),ebr,forcon(i)
 198         enddo
 199       endif
 200       write (iout,'(a)')
 201       return
 202    11 stop "Error reading reference structure"
 203       end
 204 c-----------------------------------------------------------------------------
 205       logical function seq_comp(itypea,itypeb,length)
 206       implicit none
 207       integer length,itypea(length),itypeb(length)
 208       integer i
 209       do i=1,length
 210         if (itypea(i).ne.itypeb(i)) then
 211           seq_comp=.false.
 212           return
 213         endif
 214       enddo
 215       seq_comp=.true.
 216       return
 217       end
 218 c-----------------------------------------------------------------------------
 219       subroutine read_bridge
 220 C Read information about disulfide bridges.
 221       implicit real*8 (a-h,o-z)
 222       include 'DIMENSIONS'
 223       include 'DIMENSIONS.ZSCOPT'
 224       include 'COMMON.IOUNITS'
 225       include 'COMMON.GEO'
 226       include 'COMMON.VAR'
 227       include 'COMMON.INTERACT'
 228       include 'COMMON.NAMES'
 229       include 'COMMON.CHAIN'
 230       include 'COMMON.FFIELD'
 231       include 'COMMON.SBRIDGE'
 232 C Read bridging residues.
 233       read (inp,*) ns,(iss(i),i=1,ns)
 234       print *,'ns=',ns
 235       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 236 C Check whether the specified bridging residues are cystines.
 237       do i=1,ns
 238         if (itype(iss(i)).ne.1) then
 239           write (iout,'(2a,i3,a)')
 240      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 241      &   ' can form a disulfide bridge?!!!'
 242           write (*,'(2a,i3,a)')
 243      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 244      &   ' can form a disulfide bridge?!!!'
 245          stop
 246         endif
 247       enddo
 248 C Read preformed bridges.
 249       if (ns.gt.0) then
 250       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 251       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 252       if (nss.gt.0) then
 253         nhpb=nss
 254 C Check if the residues involved in bridges are in the specified list of
 255 C bridging residues.
 256         do i=1,nss
 257           do j=1,i-1
 258             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 259      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 260               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 261      &      ' contains residues present in other pairs.'
 262               write (*,'(a,i3,a)') 'Disulfide pair',i,
 263      &      ' contains residues present in other pairs.'
 264               stop
 265             endif
 266           enddo
 267           do j=1,ns
 268             if (ihpb(i).eq.iss(j)) goto 10
 269           enddo
 270           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 271    10     continue
 272           do j=1,ns
 273             if (jhpb(i).eq.iss(j)) goto 20
 274           enddo
 275           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 276    20     continue
 277           dhpb(i)=dbr
 278           forcon(i)=fbr
 279         enddo
 280         do i=1,nss
 281           ihpb(i)=ihpb(i)+nres
 282           jhpb(i)=jhpb(i)+nres
 283         enddo
 284       endif
 285       endif
 286       if (ns.gt.0.and.dyn_ss) then
 287 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
 288 C the bond
 289           do i=nss+1,nhpb
 290 C /06/28/2013 Adasko: nss number of full SS bonds
 291             ihpb(i-nss)=ihpb(i)
 292             jhpb(i-nss)=jhpb(i)
 293             forcon(i-nss)=forcon(i)
 294             dhpb(i-nss)=dhpb(i)
 295           enddo
 296           nhpb=nhpb-nss
 297           nss=0
 298           call hpb_partition
 299           do i=1,ns
 300             dyn_ss_mask(iss(i))=.true.
 301 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
 302 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 303           enddo
 304       endif
 305       return
 306       end
 307 c------------------------------------------------------------------------------
 308       subroutine read_angles(kanal,iscor,energ,iprot,*)
 309       implicit real*8 (a-h,o-z)
 310       include 'DIMENSIONS'
 311       include 'DIMENSIONS.ZSCOPT'
 312       include 'COMMON.INTERACT'
 313       include 'COMMON.SBRIDGE'
 314       include 'COMMON.GEO'
 315       include 'COMMON.VAR'
 316       include 'COMMON.CHAIN'
 317       include 'COMMON.IOUNITS'
 318       character*80 lineh
 319       read(kanal,'(a80)',end=10,err=10) lineh
 320       read(lineh(:5),*,err=8) ic
 321       read(lineh(6:),*,err=8) energ
 322       goto 9
 323     8 ic=1
 324       print *,'error, assuming e=1d10',lineh
 325       energ=1d10
 326       nss=0
 327     9 continue
 328       read(lineh(18:),*,end=10,err=10) nss
 329       IF (NSS.LT.9) THEN
 330         read (lineh(20:),*,end=10,err=10)
 331      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 332       ELSE
 333         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 334         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 335      &    I=9,NSS),iscor
 336       ENDIF
 337 c      print *,"energy",energ," iscor",iscor
 338       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 339       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 340       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 341       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 342       do i=1,nres
 343         theta(i)=deg2rad*theta(i)
 344         phi(i)=deg2rad*phi(i)
 345         alph(i)=deg2rad*alph(i)
 346         omeg(i)=deg2rad*omeg(i)
 347       enddo
 348       return
 349    10 return1
 350       end
 351 c-------------------------------------------------------------------------------
 352       subroutine read_dist_constr
 353       implicit real*8 (a-h,o-z)
 354       include 'DIMENSIONS'
 355       include 'DIMENSIONS.ZSCOPT'
 356       include 'COMMON.CONTROL'
 357       include 'COMMON.CHAIN'
 358       include 'COMMON.IOUNITS'
 359       include 'COMMON.SBRIDGE'
 360       integer ifrag_(2,100),ipair_(2,100)
 361       double precision wfrag_(100),wpair_(100)
 362       character*500 controlcard
 363 c      write (iout,*) "Calling read_dist_constr"
 364 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 365 c      call flush(iout)
 366       call card_concat(controlcard,.true.)
 367       call readi(controlcard,"NFRAG",nfrag_,0)
 368       call readi(controlcard,"NPAIR",npair_,0)
 369       call readi(controlcard,"NDIST",ndist_,0)
 370       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 371       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 372       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 373       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 374       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 375       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 376       write (iout,*) "IFRAG"
 377       do i=1,nfrag_
 378         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 379       enddo
 380       write (iout,*) "IPAIR"
 381       do i=1,npair_
 382         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 383       enddo
 384       call flush(iout)
 385       if (.not.refstr .and. nfrag_.gt.0) then
 386         write (iout,*)
 387      &  "ERROR: no reference structure to compute distance restraints"
 388         write (iout,*)
 389      &  "Restraints must be specified explicitly (NDIST=number)"
 390         stop
 391       endif
 392       if (nfrag_.lt.2 .and. npair_.gt.0) then
 393         write (iout,*) "ERROR: Less than 2 fragments specified",
 394      &   " but distance restraints between pairs requested"
 395         stop
 396       endif
 397       call flush(iout)
 398       do i=1,nfrag_
 399         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 400         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 401      &    ifrag_(2,i)=nstart_sup+nsup-1
 402 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 403         call flush(iout)
 404         if (wfrag_(i).gt.0.0d0) then
 405         do j=ifrag_(1,i),ifrag_(2,i)-1
 406           do k=j+1,ifrag_(2,i)
 407             write (iout,*) "j",j," k",k
 408             ddjk=dist(j,k)
 409             if (constr_dist.eq.1) then
 410             nhpb=nhpb+1
 411             ihpb(nhpb)=j
 412             jhpb(nhpb)=k
 413               dhpb(nhpb)=ddjk
 414             forcon(nhpb)=wfrag_(i)
 415             else if (constr_dist.eq.2) then
 416               if (ddjk.le.dist_cut) then
 417                 nhpb=nhpb+1
 418                 ihpb(nhpb)=j
 419                 jhpb(nhpb)=k
 420                 dhpb(nhpb)=ddjk
 421                 forcon(nhpb)=wfrag_(i)
 422               endif
 423             else
 424               nhpb=nhpb+1
 425               ihpb(nhpb)=j
 426               jhpb(nhpb)=k
 427               dhpb(nhpb)=ddjk
 428               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 429             endif
 430             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 431      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 432           enddo
 433         enddo
 434         endif
 435       enddo
 436       do i=1,npair_
 437         if (wpair_(i).gt.0.0d0) then
 438         ii = ipair_(1,i)
 439         jj = ipair_(2,i)
 440         if (ii.gt.jj) then
 441           itemp=ii
 442           ii=jj
 443           jj=itemp
 444         endif
 445         do j=ifrag_(1,ii),ifrag_(2,ii)
 446           do k=ifrag_(1,jj),ifrag_(2,jj)
 447             nhpb=nhpb+1
 448             ihpb(nhpb)=j
 449             jhpb(nhpb)=k
 450             forcon(nhpb)=wpair_(i)
 451             dhpb(nhpb)=dist(j,k)
 452             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 453      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 454           enddo
 455         enddo
 456         endif
 457       enddo
 458       do i=1,ndist_
 459         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 460      &     ibecarb(i),forcon(nhpb+1)
 461         if (forcon(nhpb+1).gt.0.0d0) then
 462           nhpb=nhpb+1
 463           if (ibecarb(i).gt.0) then
 464             ihpb(i)=ihpb(i)+nres
 465             jhpb(i)=jhpb(i)+nres
 466           endif
 467           if (dhpb(nhpb).eq.0.0d0)
 468      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 469         endif
 470       enddo
 471       do i=1,nhpb
 472           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 473      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 474       enddo
 475       call flush(iout)
 476       return
 477       end
 478
 479
 480
 481 c====-------------------------------------------------------------------
 482       subroutine read_constr_homology
 483
 484       include 'DIMENSIONS'
 485       include 'DIMENSIONS.ZSCOPT'
 486 #ifdef MPI
 487       include 'mpif.h'
 488 #endif
 489       include 'COMMON.SETUP'
 490       include 'COMMON.CONTROL'
 491       include 'COMMON.CHAIN'
 492       include 'COMMON.IOUNITS'
 493       include 'COMMON.GEO'
 494       include 'COMMON.INTERACT'
 495       include 'COMMON.HOMRESTR'
 496 c
 497 c For new homol impl
 498 c
 499       include 'COMMON.VAR'
 500 c     include 'include_unres/COMMON.VAR'
 501 c
 502
 503 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 504 c    &                 dist_cut
 505 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 506 c    &    sigma_odl_temp(maxres,maxres,max_template)
 507       character*2 kic2
 508       character*24 model_ki_dist, model_ki_angle
 509       character*500 controlcard
 510       integer ki, i, j, k, l
 511       logical lprn /.true./
 512 c
 513 c     FP - Nov. 2014 Temporary specifications for new vars
 514 c
 515       double precision rescore_tmp,x12,y12,z12
 516       double precision, dimension (max_template,maxres) :: rescore
 517       character*24 tpl_k_rescore
 518 c -----------------------------------------------------------------
 519 c Reading multiple PDB ref structures and calculation of retraints
 520 c not using pre-computed ones stored in files model_ki_{dist,angle}
 521 c FP (Nov., 2014)
 522 c -----------------------------------------------------------------
 523 c
 524 c
 525 c Alternative: reading from input
 526       call card_concat(controlcard,.true.)
 527       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 528       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 529       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 530       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 531       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 532
 533       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
 534       if (homol_nset.gt.1)then
 535          call card_concat(controlcard,.true.)
 536          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 537          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 538           write(iout,*) "iset homology_weight "
 539 c         do i=1,homol_nset
 540 c          write(iout,*) i,waga_homology(i)
 541 c         enddo
 542          endif
 543          iset=mod(kolor,homol_nset)+1
 544       else
 545        iset=1
 546        waga_homology(1)=1.0
 547       endif
 548 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 549
 550 cd      write (iout,*) "nnt",nnt," nct",nct
 551 cd      call flush(iout)
 552
 553
 554       lim_odl=0
 555       lim_dih=0
 556 c
 557 c  New
 558 c
 559       lim_theta=0
 560       lim_xx=0
 561 c
 562 c  Reading HM global scores (prob not required)
 563 c
 564 c     open (4,file="HMscore")
 565 c     do k=1,constr_homology
 566 c       read (4,*,end=521) hmscore_tmp
 567 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 568 c       write(*,*) "Model", k, ":", hmscore(k)
 569 c     enddo
 570 c521  continue
 571
 572 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 573
 574       do k=1,constr_homology
 575
 576         read(inp,'(a)') pdbfile
 577 c  Next stament causes error upon compilation (?)
 578 c       if(me.eq.king.or. .not. out1file)
 579 c         write (iout,'(2a)') 'PDB data will be read from file ',
 580 c    &   pdbfile(:ilen(pdbfile))
 581         open(ipdbin,file=pdbfile,status='old',err=33)
 582         goto 34
 583   33    write (iout,'(a)') 'Error opening PDB file.'
 584         stop
 585   34    continue
 586 c        print *,'Begin reading pdb data'
 587 c
 588 c Files containing res sim or local scores (former containing sigmas)
 589 c
 590
 591         write(kic2,'(bz,i2.2)') k
 592
 593         tpl_k_rescore="template"//kic2//".sco"
 594 c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
 595 c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
 596 c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
 597 c       tpl_k_sigma_d="template"//kic2//".sigma_d"
 598
 599         unres_pdb=.false.
 600         call readpdb
 601         do i=1,2*nres
 602           do j=1,3
 603             crefjlee(j,i)=c(j,i)
 604           enddo
 605         enddo
 606 #ifdef DEBUG
 607         do i=1,nres
 608           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 609      &      (crefjlee(j,i+nres),j=1,3)
 610         enddo
 611 #endif
 612         write (iout,*) "read_constr_homology: after reading pdb file"
 613         call flush(iout)
 614
 615 c
 616 c     Distance restraints
 617 c
 618 c          ... --> odl(k,ii)
 619 C Copy the coordinates from reference coordinates (?)
 620         do i=1,2*nres
 621           do j=1,3
 622             c(j,i)=cref(j,i)
 623 c           write (iout,*) "c(",j,i,") =",c(j,i)
 624           enddo
 625         enddo
 626 c
 627 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 628 c
 629 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 630           open (ientin,file=tpl_k_rescore,status='old')
 631           do irec=1,maxdim ! loop for reading res sim
 632             if (irec.eq.1) then
 633                rescore(k,irec)=0.0d0
 634                goto 1301
 635             endif
 636             read (ientin,*,end=1401) rescore_tmp
 637 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 638             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 639 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 640  1301     continue
 641           enddo
 642  1401   continue
 643           close (ientin)
 644 c         open (ientin,file=tpl_k_sigma_odl,status='old')
 645 c         do irec=1,maxdim ! loop for reading sigma_odl
 646 c            read (ientin,*,end=1401) i, j,
 647 c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
 648 c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
 649 c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k)
 650 c         enddo
 651 c 1401   continue
 652 c         close (ientin)
 653         if (waga_dist.ne.0.0d0) then
 654           ii=0
 655           do i = nnt,nct-2 ! right? without parallel.
 656             do j=i+2,nct ! right?
 657 c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
 658 c           do j=i+2,nres ! ibid
 659 c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
 660 c           do j=i+2,nct ! ibid
 661               ii=ii+1
 662 c             write (iout,*) "k",k
 663 c             write (iout,*) "i",i," j",j," constr_homology",
 664 c    &                       constr_homology
 665               ires_homo(ii)=i
 666               jres_homo(ii)=j
 667 c
 668 c Attempt to replace dist(i,j) by its definition in ...
 669 c
 670               x12=c(1,i)-c(1,j)
 671               y12=c(2,i)-c(2,j)
 672               z12=c(3,i)-c(3,j)
 673               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 674               odl(k,ii)=distal
 675 c
 676 c             odl(k,ii)=dist(i,j)
 677 c             write (iout,*) "dist(",i,j,") =",dist(i,j)
 678 c             write (iout,*) "distal = ",distal
 679 c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
 680 c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 681 c    &                      "rescore(",k,j,") =",rescore(k,j)
 682 c
 683 c  Calculation of sigma from res sim
 684 c
 685 c             if (odl(k,ii).le.6.0d0) then
 686 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
 687 c    Other functional forms possible depending on odl(k,ii), eg.
 688 c
 689             if (odl(k,ii).le.dist_cut) then
 690               sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
 691 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
 692             else
 693               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
 694      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 695
 696 c   Following expr replaced by a positive exp argument
 697 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 698 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
 699
 700 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
 701 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
 702             endif
 703 c
 704               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
 705 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
 706 c
 707 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
 708 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
 709             enddo
 710 c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
 711           enddo
 712 c         lim_odl=ii
 713 c         if (constr_homology.gt.0) call homology_partition
 714         endif
 715 c
 716 c     Theta, dihedral and SC retraints
 717 c
 718         if (waga_angle.gt.0.0d0) then
 719 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 720 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 721 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 722 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 723 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 724 c    &                            sigma_dih(k,i+nnt-1)
 725 c         enddo
 726 c1402   continue
 727 c         close (ientin)
 728           do i = nnt+3,nct ! right? without parallel.
 729 c         do i=1,nres ! alternative for bounds acc to readpdb?
 730 c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
 731 c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
 732             dih(k,i)=phiref(i) ! right?
 733 c           read (ientin,*) sigma_dih(k,i) ! original variant
 734 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 735 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 736 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 737 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 738 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 739
 740             sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
 741      &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
 742 c
 743 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 744 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 745 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 746 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 747 c   Instead of res sim other local measure of b/b str reliability possible
 748             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 749 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 750             if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
 751 c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
 752           enddo
 753         endif
 754
 755         if (waga_theta.gt.0.0d0) then
 756 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 757 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 758 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 759 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 760 c    &                              sigma_theta(k,i+nnt-1)
 761 c         enddo
 762 c1403   continue
 763 c         close (ientin)
 764
 765           do i = nnt+2,nct ! right? without parallel.
 766 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 767 c         do i=ithet_start,ithet_end ! with FG parallel.
 768              thetatpl(k,i)=thetaref(i)
 769 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 770 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 771 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 772 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 773 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 774              sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
 775      &                        rescore(k,i-2) !  right expression ?
 776              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 777
 778 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 779 c                             rescore(k,i-2) !  right expression ?
 780 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 781              if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
 782           enddo
 783         endif
 784
 785         if (waga_d.gt.0.0d0) then
 786 c       open (ientin,file=tpl_k_sigma_d,status='old')
 787 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 788 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 789 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 790 c    &                          sigma_d(k,i+nnt-1)
 791 c         enddo
 792 c1404   continue
 793           close (ientin)
 794
 795           do i = nnt,nct ! right? without parallel.
 796 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 797 c         do i=loc_start,loc_end ! with FG parallel.
 798              if (itype(i).eq.10) goto 1 ! right?
 799                xxtpl(k,i)=xxref(i)
 800                yytpl(k,i)=yyref(i)
 801                zztpl(k,i)=zzref(i)
 802 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 803 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 804 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 805 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 806                sigma_d(k,i)=rescore(k,i) !  right expression ?
 807                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 808
 809 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 810 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 811 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 812                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
 813     1     continue
 814           enddo
 815         endif
 816         close(ientin)
 817       enddo
 818       if (waga_dist.ne.0.0d0) lim_odl=ii
 819       if (constr_homology.gt.0) call homology_partition
 820       if (constr_homology.gt.0) call init_int_table
 821 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 822 cd     & "lim_xx=",lim_xx
 823 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 824 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 825 c
 826 c Print restraints
 827 c
 828       if (.not.lprn) return
 829 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 830       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 831        write (iout,*) "Distance restraints from templates"
 832        do ii=1,lim_odl
 833        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
 834      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
 835        enddo
 836        write (iout,*) "Dihedral angle restraints from templates"
 837        do i=nnt+3,lim_dih
 838         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
 839      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 840        enddo
 841        write (iout,*) "Virtual-bond angle restraints from templates"
 842        do i=nnt+2,lim_theta
 843         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
 844      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 845        enddo
 846        write (iout,*) "SC restraints from templates"
 847        do i=nnt,lim_xx
 848         write(iout,'(i5,10(4f8.2,4x))') i,
 849      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 850      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 851        enddo
 852       endif
 853 c -----------------------------------------------------------------
 854       return
 855       end
 856 c----------------------------------------------------------------------