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