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