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 #ifdef OLDSIGMA
 694               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
 695      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 696 #else
 697               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
 698      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 699 #endif
 700 c   Following expr replaced by a positive exp argument
 701 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 702 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
 703
 704 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
 705 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
 706             endif
 707 c
 708               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
 709 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
 710 c
 711 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
 712 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
 713             enddo
 714 c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
 715           enddo
 716 c         lim_odl=ii
 717 c         if (constr_homology.gt.0) call homology_partition
 718         endif
 719 c
 720 c     Theta, dihedral and SC retraints
 721 c
 722         if (waga_angle.gt.0.0d0) then
 723 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 724 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 725 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 726 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 727 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 728 c    &                            sigma_dih(k,i+nnt-1)
 729 c         enddo
 730 c1402   continue
 731 c         close (ientin)
 732           do i = nnt+3,nct ! right? without parallel.
 733 c         do i=1,nres ! alternative for bounds acc to readpdb?
 734 c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
 735 c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
 736             dih(k,i)=phiref(i) ! right?
 737 c           read (ientin,*) sigma_dih(k,i) ! original variant
 738 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 739 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 740 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 741 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 742 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 743
 744             sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
 745      &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
 746 c
 747 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 748 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 749 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 750 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 751 c   Instead of res sim other local measure of b/b str reliability possible
 752             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 753 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 754             if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
 755 c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
 756           enddo
 757         endif
 758
 759         if (waga_theta.gt.0.0d0) then
 760 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 761 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 762 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 763 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 764 c    &                              sigma_theta(k,i+nnt-1)
 765 c         enddo
 766 c1403   continue
 767 c         close (ientin)
 768
 769           do i = nnt+2,nct ! right? without parallel.
 770 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 771 c         do i=ithet_start,ithet_end ! with FG parallel.
 772              thetatpl(k,i)=thetaref(i)
 773 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 774 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 775 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 776 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 777 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 778              sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
 779      &                        rescore(k,i-2) !  right expression ?
 780              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 781
 782 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 783 c                             rescore(k,i-2) !  right expression ?
 784 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 785              if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
 786           enddo
 787         endif
 788
 789         if (waga_d.gt.0.0d0) then
 790 c       open (ientin,file=tpl_k_sigma_d,status='old')
 791 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 792 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 793 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 794 c    &                          sigma_d(k,i+nnt-1)
 795 c         enddo
 796 c1404   continue
 797           close (ientin)
 798
 799           do i = nnt,nct ! right? without parallel.
 800 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 801 c         do i=loc_start,loc_end ! with FG parallel.
 802              if (itype(i).eq.10) goto 1 ! right?
 803                xxtpl(k,i)=xxref(i)
 804                yytpl(k,i)=yyref(i)
 805                zztpl(k,i)=zzref(i)
 806 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 807 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 808 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 809 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 810                sigma_d(k,i)=rescore(k,i) !  right expression ?
 811                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 812
 813 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 814 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 815 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 816                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
 817     1     continue
 818           enddo
 819         endif
 820         close(ientin)
 821       enddo
 822       if (waga_dist.ne.0.0d0) lim_odl=ii
 823       if (constr_homology.gt.0) call homology_partition
 824       if (constr_homology.gt.0) call init_int_table
 825 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 826 cd     & "lim_xx=",lim_xx
 827 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 828 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 829 c
 830 c Print restraints
 831 c
 832       if (.not.lprn) return
 833 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 834       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 835        write (iout,*) "Distance restraints from templates"
 836        do ii=1,lim_odl
 837        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
 838      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
 839        enddo
 840        write (iout,*) "Dihedral angle restraints from templates"
 841        do i=nnt+3,lim_dih
 842         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
 843      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 844        enddo
 845        write (iout,*) "Virtual-bond angle restraints from templates"
 846        do i=nnt+2,lim_theta
 847         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
 848      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 849        enddo
 850        write (iout,*) "SC restraints from templates"
 851        do i=nnt,lim_xx
 852         write(iout,'(i5,10(4f8.2,4x))') i,
 853      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 854      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 855        enddo
 856       endif
 857 c -----------------------------------------------------------------
 858       return
 859       end
 860 c----------------------------------------------------------------------