source/wham/src-old/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       logical unres_pdb
 513 c
 514 c     FP - Nov. 2014 Temporary specifications for new vars
 515 c
 516       double precision rescore_tmp,x12,y12,z12
 517       double precision, dimension (max_template,maxres) :: rescore
 518       character*24 tpl_k_rescore
 519 c -----------------------------------------------------------------
 520 c Reading multiple PDB ref structures and calculation of retraints
 521 c not using pre-computed ones stored in files model_ki_{dist,angle}
 522 c FP (Nov., 2014)
 523 c -----------------------------------------------------------------
 524 c
 525 c
 526 c Alternative: reading from input
 527       call card_concat(controlcard,.true.)
 528       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 529       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 530       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 531       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 532       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 533
 534       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
 535       if (homol_nset.gt.1)then
 536          call card_concat(controlcard,.true.)
 537          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 538          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 539           write(iout,*) "iset homology_weight "
 540 c         do i=1,homol_nset
 541 c          write(iout,*) i,waga_homology(i)
 542 c         enddo
 543          endif
 544          iset=mod(kolor,homol_nset)+1
 545       else
 546        iset=1
 547        waga_homology(1)=1.0
 548       endif
 549 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 550
 551 cd      write (iout,*) "nnt",nnt," nct",nct
 552 cd      call flush(iout)
 553
 554
 555       lim_odl=0
 556       lim_dih=0
 557 c
 558 c  New
 559 c
 560       lim_theta=0
 561       lim_xx=0
 562 c
 563 c  Reading HM global scores (prob not required)
 564 c
 565 c     open (4,file="HMscore")
 566 c     do k=1,constr_homology
 567 c       read (4,*,end=521) hmscore_tmp
 568 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 569 c       write(*,*) "Model", k, ":", hmscore(k)
 570 c     enddo
 571 c521  continue
 572
 573 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 574
 575       do k=1,constr_homology
 576
 577         read(inp,'(a)') pdbfile
 578 c  Next stament causes error upon compilation (?)
 579 c       if(me.eq.king.or. .not. out1file)
 580 c         write (iout,'(2a)') 'PDB data will be read from file ',
 581 c    &   pdbfile(:ilen(pdbfile))
 582         open(ipdbin,file=pdbfile,status='old',err=33)
 583         goto 34
 584   33    write (iout,'(a)') 'Error opening PDB file.'
 585         stop
 586   34    continue
 587 c        print *,'Begin reading pdb data'
 588 c
 589 c Files containing res sim or local scores (former containing sigmas)
 590 c
 591
 592         write(kic2,'(bz,i2.2)') k
 593
 594         tpl_k_rescore="template"//kic2//".sco"
 595 c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
 596 c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
 597 c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
 598 c       tpl_k_sigma_d="template"//kic2//".sigma_d"
 599
 600         unres_pdb=.false.
 601         call readpdb
 602         do i=1,2*nres
 603           do j=1,3
 604             crefjlee(j,i)=c(j,i)
 605           enddo
 606         enddo
 607 #ifdef DEBUG
 608         do i=1,nres
 609           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 610      &      (crefjlee(j,i+nres),j=1,3)
 611         enddo
 612 #endif
 613         write (iout,*) "read_constr_homology: after reading pdb file"
 614         call flush(iout)
 615
 616 c
 617 c     Distance restraints
 618 c
 619 c          ... --> odl(k,ii)
 620 C Copy the coordinates from reference coordinates (?)
 621         do i=1,2*nres
 622           do j=1,3
 623             c(j,i)=cref(j,i)
 624 c           write (iout,*) "c(",j,i,") =",c(j,i)
 625           enddo
 626         enddo
 627 c
 628 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 629 c
 630 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 631           open (ientin,file=tpl_k_rescore,status='old')
 632           do irec=1,maxdim ! loop for reading res sim
 633             if (irec.eq.1) then
 634                rescore(k,irec)=0.0d0
 635                goto 1301
 636             endif
 637             read (ientin,*,end=1401) rescore_tmp
 638 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 639             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 640 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 641  1301     continue
 642           enddo
 643  1401   continue
 644           close (ientin)
 645 c         open (ientin,file=tpl_k_sigma_odl,status='old')
 646 c         do irec=1,maxdim ! loop for reading sigma_odl
 647 c            read (ientin,*,end=1401) i, j,
 648 c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
 649 c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
 650 c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k)
 651 c         enddo
 652 c 1401   continue
 653 c         close (ientin)
 654         if (waga_dist.ne.0.0d0) then
 655           ii=0
 656           do i = nnt,nct-2 ! right? without parallel.
 657             do j=i+2,nct ! right?
 658 c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
 659 c           do j=i+2,nres ! ibid
 660 c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
 661 c           do j=i+2,nct ! ibid
 662               ii=ii+1
 663 c             write (iout,*) "k",k
 664 c             write (iout,*) "i",i," j",j," constr_homology",
 665 c    &                       constr_homology
 666               ires_homo(ii)=i
 667               jres_homo(ii)=j
 668 c
 669 c Attempt to replace dist(i,j) by its definition in ...
 670 c
 671               x12=c(1,i)-c(1,j)
 672               y12=c(2,i)-c(2,j)
 673               z12=c(3,i)-c(3,j)
 674               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 675               odl(k,ii)=distal
 676 c
 677 c             odl(k,ii)=dist(i,j)
 678 c             write (iout,*) "dist(",i,j,") =",dist(i,j)
 679 c             write (iout,*) "distal = ",distal
 680 c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
 681 c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 682 c    &                      "rescore(",k,j,") =",rescore(k,j)
 683 c
 684 c  Calculation of sigma from res sim
 685 c
 686 c             if (odl(k,ii).le.6.0d0) then
 687 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
 688 c    Other functional forms possible depending on odl(k,ii), eg.
 689 c
 690             if (odl(k,ii).le.dist_cut) then
 691               sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
 692 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
 693             else
 694 #ifdef OLDSIGMA
 695               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
 696      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 697 #else
 698               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
 699      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 700 #endif
 701 c   Following expr replaced by a positive exp argument
 702 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 703 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
 704
 705 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
 706 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
 707             endif
 708 c
 709               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
 710 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
 711 c
 712 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
 713 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
 714             enddo
 715 c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
 716           enddo
 717 c         lim_odl=ii
 718 c         if (constr_homology.gt.0) call homology_partition
 719         endif
 720 c
 721 c     Theta, dihedral and SC retraints
 722 c
 723         if (waga_angle.gt.0.0d0) then
 724 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 725 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 726 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 727 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 728 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 729 c    &                            sigma_dih(k,i+nnt-1)
 730 c         enddo
 731 c1402   continue
 732 c         close (ientin)
 733           do i = nnt+3,nct ! right? without parallel.
 734 c         do i=1,nres ! alternative for bounds acc to readpdb?
 735 c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
 736 c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
 737             dih(k,i)=phiref(i) ! right?
 738 c           read (ientin,*) sigma_dih(k,i) ! original variant
 739 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 740 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 741 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 742 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 743 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 744
 745             sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
 746      &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
 747 c
 748 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 749 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 750 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 751 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 752 c   Instead of res sim other local measure of b/b str reliability possible
 753             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 754 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 755             if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
 756 c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
 757           enddo
 758         endif
 759
 760         if (waga_theta.gt.0.0d0) then
 761 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 762 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 763 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 764 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 765 c    &                              sigma_theta(k,i+nnt-1)
 766 c         enddo
 767 c1403   continue
 768 c         close (ientin)
 769
 770           do i = nnt+2,nct ! right? without parallel.
 771 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 772 c         do i=ithet_start,ithet_end ! with FG parallel.
 773              thetatpl(k,i)=thetaref(i)
 774 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 775 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 776 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 777 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 778 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 779              sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
 780      &                        rescore(k,i-2) !  right expression ?
 781              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 782
 783 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 784 c                             rescore(k,i-2) !  right expression ?
 785 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 786              if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
 787           enddo
 788         endif
 789
 790         if (waga_d.gt.0.0d0) then
 791 c       open (ientin,file=tpl_k_sigma_d,status='old')
 792 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 793 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 794 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 795 c    &                          sigma_d(k,i+nnt-1)
 796 c         enddo
 797 c1404   continue
 798           close (ientin)
 799
 800           do i = nnt,nct ! right? without parallel.
 801 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 802 c         do i=loc_start,loc_end ! with FG parallel.
 803              if (itype(i).eq.10) goto 1 ! right?
 804                xxtpl(k,i)=xxref(i)
 805                yytpl(k,i)=yyref(i)
 806                zztpl(k,i)=zzref(i)
 807 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 808 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 809 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 810 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 811                sigma_d(k,i)=rescore(k,i) !  right expression ?
 812                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 813
 814 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 815 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 816 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 817                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
 818     1     continue
 819           enddo
 820         endif
 821         close(ientin)
 822       enddo
 823       if (waga_dist.ne.0.0d0) lim_odl=ii
 824       if (constr_homology.gt.0) call homology_partition
 825       if (constr_homology.gt.0) call init_int_table
 826 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 827 cd     & "lim_xx=",lim_xx
 828 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 829 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 830 c
 831 c Print restraints
 832 c
 833       if (.not.lprn) return
 834 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 835       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 836        write (iout,*) "Distance restraints from templates"
 837        do ii=1,lim_odl
 838        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
 839      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
 840        enddo
 841        write (iout,*) "Dihedral angle restraints from templates"
 842        do i=nnt+3,lim_dih
 843         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
 844      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 845        enddo
 846        write (iout,*) "Virtual-bond angle restraints from templates"
 847        do i=nnt+2,lim_theta
 848         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
 849      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 850        enddo
 851        write (iout,*) "SC restraints from templates"
 852        do i=nnt,lim_xx
 853         write(iout,'(i5,10(4f8.2,4x))') i,
 854      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 855      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 856        enddo
 857       endif
 858 c -----------------------------------------------------------------
 859       return
 860       end
 861 c----------------------------------------------------------------------