source/cluster/wham/src/readrtns.F

   1       subroutine read_control
   2 C
   3 C Read molecular data
   4 C
   5       implicit none
   6       include 'DIMENSIONS'
   7       include 'sizesclu.dat'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.TIME1'
  10       include 'COMMON.SBRIDGE'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.HEADER'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.FREE'
  17       character*320 controlcard,ucase
  18 #ifdef MPL
  19       include 'COMMON.INFO'
  20 #endif
  21       integer i
  22
  23       read (INP,'(a80)') titel
  24       call card_concat(controlcard)
  25
  26       call readi(controlcard,'NRES',nres,0)
  27       call readi(controlcard,'RESCALE',rescale_mode,2)
  28       call readi(controlcard,'PDBOUT',outpdb,0)
  29       call readi(controlcard,'MOL2OUT',outmol2,0)
  30       refstr=(index(controlcard,'REFSTR').gt.0)
  31       write (iout,*) "REFSTR",refstr
  32       pdbref=(index(controlcard,'PDBREF').gt.0)
  33       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
  34       iscode=index(controlcard,'ONE_LETTER')
  35       tree=(index(controlcard,'MAKE_TREE').gt.0)
  36       min_var=(index(controlcard,'MINVAR').gt.0)
  37       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
  38       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
  39       call readi(controlcard,'NCUT',ncut,0)
  40       call readi(controlcard,'NCLUST',nclust,5)
  41       call readi(controlcard,'NSTART',nstart,0)
  42       call readi(controlcard,'NEND',nend,0)
  43       call reada(controlcard,'ECUT',ecut,10.0d0)
  44       call reada(controlcard,'PROB',prob_limit,0.99d0)
  45       write (iout,*) "Probability limit",prob_limit
  46       lgrp=(index(controlcard,'LGRP').gt.0)
  47       caonly=(index(controlcard,'CA_ONLY').gt.0)
  48       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
  49       if (ncut.gt.0)
  50      & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
  51       call readi(controlcard,'IOPT',iopt,2)
  52       lside = index(controlcard,"SIDE").gt.0
  53       efree = index(controlcard,"EFREE").gt.0
  54       call readi(controlcard,'NTEMP',nT,1)
  55       write (iout,*) "nT",nT
  56       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
  57       write (iout,*) "nT",nT
  58       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
  59       do i=1,nT
  60         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
  61       enddo
  62       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
  63       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
  64       lprint_int=index(controlcard,"PRINT_INT") .gt.0
  65       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  66       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  67       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  68      & " CONSTR_DIST",constr_dist
  69 <<<<<<< HEAD
  70 =======
  71
  72       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
  73       write (iout,*) "with_homology_constr ",with_dihed_constr,
  74      & " CONSTR_HOMOLOGY",constr_homology
  75       print_homology_restraints=
  76      & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
  77       print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
  78       print_homology_models=
  79      & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
  80
  81 #ifdef AIX
  82       call flush_(iout)
  83 #else
  84 >>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit
  85       call flush(iout)
  86       if (min_var) iopt=1
  87       return
  88       end
  89 c--------------------------------------------------------------------------
  90       subroutine molread
  91 C
  92 C Read molecular data.
  93 C
  94       implicit none
  95       include 'DIMENSIONS'
  96       include 'COMMON.IOUNITS'
  97       include 'COMMON.GEO'
  98       include 'COMMON.VAR'
  99       include 'COMMON.INTERACT'
 100       include 'COMMON.LOCAL'
 101       include 'COMMON.NAMES'
 102       include 'COMMON.CHAIN'
 103       include 'COMMON.FFIELD'
 104       include 'COMMON.SBRIDGE'
 105       include 'COMMON.HEADER'
 106       include 'COMMON.CONTROL'
 107       include 'COMMON.CONTACTS'
 108       include 'COMMON.TIME1'
 109       include 'COMMON.TORCNSTR'
 110 #ifdef MPL
 111       include 'COMMON.INFO'
 112 #endif
 113       character*4 sequence(maxres)
 114       character*800 weightcard
 115       integer rescode
 116       double precision x(maxvar)
 117       integer itype_pdb(maxres)
 118       logical seq_comp
 119       integer i,j
 120 C
 121 C Body
 122 C
 123 C Read weights of the subsequent energy terms.
 124       call card_concat(weightcard)
 125       call reada(weightcard,'WSC',wsc,1.0d0)
 126       call reada(weightcard,'WLONG',wsc,wsc)
 127       call reada(weightcard,'WSCP',wscp,1.0d0)
 128       call reada(weightcard,'WELEC',welec,1.0D0)
 129       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 130       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 131       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 132       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 133       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 134       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 135       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 136       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 137       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 138       call reada(weightcard,'WBOND',wbond,1.0D0)
 139       call reada(weightcard,'WTOR',wtor,1.0D0)
 140       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 141       call reada(weightcard,'WANG',wang,1.0D0)
 142       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 143       call reada(weightcard,'SCAL14',scal14,0.4D0)
 144       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 145       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 146       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 147       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 148       call reada(weightcard,"D0CM",d0cm,3.78d0)
 149       call reada(weightcard,"AKCM",akcm,15.1d0)
 150       call reada(weightcard,"AKTH",akth,11.0d0)
 151       call reada(weightcard,"AKCT",akct,12.0d0)
 152       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 153       call reada(weightcard,"V2SS",v2ss,7.61d0)
 154       call reada(weightcard,"V3SS",v3ss,13.7d0)
 155       call reada(weightcard,"EBR",ebr,-5.50D0)
 156       if (index(weightcard,'SOFT').gt.0) ipot=6
 157 C 12/1/95 Added weight for the multi-body term WCORR
 158       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 159       do i=1,maxres-1
 160         do j=i+1,maxres
 161           dyn_ssbond_ij(i,j)=1.0d300
 162         enddo
 163       enddo
 164       call reada(weightcard,"HT",Ht,0.0D0)
 165       if (dyn_ss) then
 166         ss_depth=ebr/wsc-0.25*eps(1,1)
 167         Ht=Ht/wsc-0.25*eps(1,1)
 168         akcm=akcm*wstrain/wsc
 169         akth=akth*wstrain/wsc
 170         akct=akct*wstrain/wsc
 171         v1ss=v1ss*wstrain/wsc
 172         v2ss=v2ss*wstrain/wsc
 173         v3ss=v3ss*wstrain/wsc
 174       else
 175         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 176       endif
 177       write (iout,'(/a)') "Disulfide bridge parameters:"
 178       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 179         write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
 180       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 181       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 182       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 183      & ' v3ss:',v3ss
 184       write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
 185       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 186       weights(1)=wsc
 187       weights(2)=wscp
 188       weights(3)=welec
 189       weights(4)=wcorr
 190       weights(5)=wcorr5
 191       weights(6)=wcorr6
 192       weights(7)=wel_loc
 193       weights(8)=wturn3
 194       weights(9)=wturn4
 195       weights(10)=wturn6
 196       weights(11)=wang
 197       weights(12)=wscloc
 198       weights(13)=wtor
 199       weights(14)=wtor_d
 200       weights(15)=wstrain
 201       weights(16)=wvdwpp
 202       weights(17)=wbond
 203       weights(18)=scal14
 204       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 205      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 206      &  wturn4,wturn6,wsccor
 207    10 format (/'Energy-term weights (unscaled):'//
 208      & 'WSCC=   ',f10.6,' (SC-SC)'/
 209      & 'WSCP=   ',f10.6,' (SC-p)'/
 210      & 'WELEC=  ',f10.6,' (p-p electr)'/
 211      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 212      & 'WBOND=  ',f10.6,' (stretching)'/
 213      & 'WANG=   ',f10.6,' (bending)'/
 214      & 'WSCLOC= ',f10.6,' (SC local)'/
 215      & 'WTOR=   ',f10.6,' (torsional)'/
 216      & 'WTORD=  ',f10.6,' (double torsional)'/
 217      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 218      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 219      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 220      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 221      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 222      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 223      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 224      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 225      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 226       if (wcorr4.gt.0.0d0) then
 227         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 228      &   'between contact pairs of peptide groups'
 229         write (iout,'(2(a,f5.3/))')
 230      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 231      &  'Range of quenching the correlation terms:',2*delt_corr
 232       else if (wcorr.gt.0.0d0) then
 233         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 234      &   'between contact pairs of peptide groups'
 235       endif
 236       write (iout,'(a,f8.3)')
 237      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 238       write (iout,'(a,f8.3)')
 239      &  'General scaling factor of SC-p interactions:',scalscp
 240       r0_corr=cutoff_corr-delt_corr
 241       do i=1,20
 242         aad(i,1)=scalscp*aad(i,1)
 243         aad(i,2)=scalscp*aad(i,2)
 244         bad(i,1)=scalscp*bad(i,1)
 245         bad(i,2)=scalscp*bad(i,2)
 246       enddo
 247
 248       call flush(iout)
 249       print *,'indpdb=',indpdb,' pdbref=',pdbref
 250
 251 C Read sequence if not taken from the pdb file.
 252       if (iscode.gt.0) then
 253         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 254       else
 255         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 256       endif
 257 C Convert sequence to numeric code
 258       do i=1,nres
 259         itype(i)=rescode(i,sequence(i),iscode)
 260       enddo
 261       if (itype(2).eq.10) then
 262         write (iout,*)
 263      &   "Glycine is the first full residue, initial dummy deleted"
 264         do i=1,nres
 265           itype(i)=itype(i+1)
 266         enddo
 267         nres=nres-1
 268       endif
 269       if (itype(nres).eq.10) then
 270         write (iout,*)
 271      &   "Glycine is the last full residue, terminal dummy deleted"
 272         nres=nres-1
 273       endif
 274       print *,nres
 275       print '(20i4)',(itype(i),i=1,nres)
 276
 277       do i=1,nres
 278 #ifdef PROCOR
 279         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 280 #else
 281         if (itype(i).eq.21) then
 282 #endif
 283           itel(i)=0
 284 #ifdef PROCOR
 285         else if (itype(i+1).ne.20) then
 286 #else
 287         else if (itype(i).ne.20) then
 288 #endif
 289           itel(i)=1
 290         else
 291           itel(i)=2
 292         endif
 293       enddo
 294       write (iout,*) "ITEL"
 295       do i=1,nres-1
 296         write (iout,*) i,itype(i),itel(i)
 297       enddo
 298
 299       print *,'Call Read_Bridge.'
 300       call read_bridge
 301       if (with_dihed_constr) then
 302
 303       read (inp,*) ndih_constr
 304       if (ndih_constr.gt.0) then
 305         read (inp,*) ftors
 306         write (iout,*) 'FTORS',ftors
 307         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 308         write (iout,*)
 309      &   'There are',ndih_constr,' constraints on phi angles.'
 310         do i=1,ndih_constr
 311           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 312         enddo
 313         do i=1,ndih_constr
 314           phi0(i)=deg2rad*phi0(i)
 315           drange(i)=deg2rad*drange(i)
 316         enddo
 317       endif
 318
 319       endif
 320       nnt=1
 321       nct=nres
 322       print *,'NNT=',NNT,' NCT=',NCT
 323       if (itype(1).eq.21) nnt=2
 324       if (itype(nres).eq.21) nct=nct-1
 325       if (nstart.lt.nnt) nstart=nnt
 326       if (nend.gt.nct .or. nend.eq.0) nend=nct
 327       write (iout,*) "nstart",nstart," nend",nend
 328       nres0=nres
 329 <<<<<<< HEAD
 330 =======
 331
 332 C     Juyong:READ init_vars
 333 C     Initialize variables!
 334 C     Juyong:READ read_info
 335 C     READ fragment information!!
 336 C     both routines should be in dfa.F file!!
 337
 338       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 339      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 340 #ifdef DEBUG
 341        write (iout,*) "Calling init_dfa_vars"
 342 #ifdef AIX
 343       call flush_(iout)
 344 #else
 345       call flush(iout)
 346 #endif
 347 #endif
 348        call init_dfa_vars
 349 #ifdef DEBUG
 350        write (iout,*) 'init_dfa_vars finished!'
 351 #ifdef AIX
 352       call flush_(iout)
 353 #else
 354       call flush(iout)
 355 #endif
 356 #endif
 357        call read_dfa_info
 358 #ifdef DEBUG
 359        write (iout,*) 'read_dfa_info finished!'
 360 #ifdef AIX
 361       call flush_(iout)
 362 #else
 363       call flush(iout)
 364 #endif
 365 #endif
 366       endif
 367
 368       if (constr_homology.gt.0) then
 369         call read_constr_homology(print_homology_restraints)
 370       endif
 371
 372 >>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit
 373 c      if (pdbref) then
 374 c        read(inp,'(a)') pdbfile
 375 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 376 c        open(ipdbin,file=pdbfile,status='old',err=33)
 377 c        goto 34
 378 c  33    write (iout,'(a)') 'Error opening PDB file.'
 379 c        stop
 380 c  34    continue
 381 c        print *,'Begin reading pdb data'
 382 c        call readpdb
 383 c        print *,'Finished reading pdb data'
 384 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 385 c        do i=1,nres
 386 c          itype_pdb(i)=itype(i)
 387 c        enddo
 388 c        close (ipdbin)
 389 c        write (iout,'(a,i3)') 'nsup=',nsup
 390 c        nstart_seq=nnt
 391 c        if (nsup.le.(nct-nnt+1)) then
 392 c          do i=0,nct-nnt+1-nsup
 393 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 394 c              nstart_seq=nnt+i
 395 c              goto 111
 396 c            endif
 397 c          enddo
 398 c          write (iout,'(a)')
 399 c     &            'Error - sequences to be superposed do not match.'
 400 c          stop
 401 c        else
 402 c          do i=0,nsup-(nct-nnt+1)
 403 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 404 c     &      then
 405 c              nstart_sup=nstart_sup+i
 406 c              nsup=nct-nnt+1
 407 c              goto 111
 408 c            endif
 409 c          enddo
 410 c          write (iout,'(a)')
 411 c     &            'Error - sequences to be superposed do not match.'
 412 c        endif
 413 c  111   continue
 414 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 415 c     &                 ' nstart_seq=',nstart_seq
 416 c      endif
 417       call init_int_table
 418       call setup_var
 419       write (iout,*) "molread: REFSTR",refstr
 420       if (refstr) then
 421         if (.not.pdbref) then
 422           call read_angles(inp,*38)
 423           goto 39
 424    38     write (iout,'(a)') 'Error reading reference structure.'
 425 #ifdef MPL
 426           call mp_stopall(Error_Msg)
 427 #else
 428           stop 'Error reading reference structure'
 429 #endif
 430    39     call chainbuild
 431           nstart_sup=nnt
 432           nstart_seq=nnt
 433           nsup=nct-nnt+1
 434           do i=1,2*nres
 435             do j=1,3
 436               cref(j,i)=c(j,i)
 437             enddo
 438           enddo
 439         endif
 440         call contact(print_contact_map,ncont_ref,icont_ref)
 441       endif
 442 c Read distance restraints
 443       if (constr_dist.gt.0) then
 444         call read_dist_constr
 445         call hpb_partition
 446       endif
 447       return
 448       end
 449 c-----------------------------------------------------------------------------
 450       logical function seq_comp(itypea,itypeb,length)
 451       implicit none
 452       integer length,itypea(length),itypeb(length)
 453       integer i
 454       do i=1,length
 455         if (itypea(i).ne.itypeb(i)) then
 456           seq_comp=.false.
 457           return
 458         endif
 459       enddo
 460       seq_comp=.true.
 461       return
 462       end
 463 c-----------------------------------------------------------------------------
 464       subroutine read_bridge
 465 C Read information about disulfide bridges.
 466       implicit none
 467       include 'DIMENSIONS'
 468       include 'COMMON.IOUNITS'
 469       include 'COMMON.GEO'
 470       include 'COMMON.VAR'
 471       include 'COMMON.INTERACT'
 472       include 'COMMON.LOCAL'
 473       include 'COMMON.NAMES'
 474       include 'COMMON.CHAIN'
 475       include 'COMMON.FFIELD'
 476       include 'COMMON.SBRIDGE'
 477       include 'COMMON.HEADER'
 478       include 'COMMON.CONTROL'
 479       include 'COMMON.TIME1'
 480 #ifdef MPL
 481       include 'COMMON.INFO'
 482 #endif
 483       integer i,j
 484 C Read bridging residues.
 485       read (inp,*) ns,(iss(i),i=1,ns)
 486       print *,'ns=',ns
 487 C Check whether the specified bridging residues are cystines.
 488       do i=1,ns
 489         if (itype(iss(i)).ne.1) then
 490           write (iout,'(2a,i3,a)')
 491      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 492      &   ' can form a disulfide bridge?!!!'
 493           write (*,'(2a,i3,a)')
 494      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 495      &   ' can form a disulfide bridge?!!!'
 496 #ifdef MPL
 497          call mp_stopall(error_msg)
 498 #else
 499          stop
 500 #endif
 501         endif
 502       enddo
 503 C Read preformed bridges.
 504       if (ns.gt.0) then
 505       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 506       if (nss.gt.0) then
 507         nhpb=nss
 508 C Check if the residues involved in bridges are in the specified list of
 509 C bridging residues.
 510         do i=1,nss
 511           do j=1,i-1
 512             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 513      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 514               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 515      &      ' contains residues present in other pairs.'
 516               write (*,'(a,i3,a)') 'Disulfide pair',i,
 517      &      ' contains residues present in other pairs.'
 518 #ifdef MPL
 519               call mp_stopall(error_msg)
 520 #else
 521               stop
 522 #endif
 523             endif
 524           enddo
 525           do j=1,ns
 526             if (ihpb(i).eq.iss(j)) goto 10
 527           enddo
 528           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 529    10     continue
 530           do j=1,ns
 531             if (jhpb(i).eq.iss(j)) goto 20
 532           enddo
 533           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 534    20     continue
 535 c          dhpb(i)=dbr
 536 c          forcon(i)=fbr
 537         enddo
 538         do i=1,nss
 539           ihpb(i)=ihpb(i)+nres
 540           jhpb(i)=jhpb(i)+nres
 541         enddo
 542       endif
 543       endif
 544       if (ns.gt.0.and.dyn_ss) then
 545           do i=nss+1,nhpb
 546             ihpb(i-nss)=ihpb(i)
 547             jhpb(i-nss)=jhpb(i)
 548             forcon(i-nss)=forcon(i)
 549             dhpb(i-nss)=dhpb(i)
 550           enddo
 551           nhpb=nhpb-nss
 552           nss=0
 553           call hpb_partition
 554           do i=1,ns
 555             dyn_ss_mask(iss(i))=.true.
 556 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 557           enddo
 558       endif
 559       print *, "Leaving brigde read"
 560       return
 561       end
 562 c----------------------------------------------------------------------------
 563       subroutine read_angles(kanal,*)
 564       implicit none
 565       include 'DIMENSIONS'
 566       include 'COMMON.GEO'
 567       include 'COMMON.VAR'
 568       include 'COMMON.CHAIN'
 569       include 'COMMON.IOUNITS'
 570       integer i,kanal
 571       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 572       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 573       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 574       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 575       do i=1,nres
 576         theta(i)=deg2rad*theta(i)
 577         phi(i)=deg2rad*phi(i)
 578         alph(i)=deg2rad*alph(i)
 579         omeg(i)=deg2rad*omeg(i)
 580       enddo
 581       return
 582    10 return1
 583       end
 584 c----------------------------------------------------------------------------
 585       subroutine reada(rekord,lancuch,wartosc,default)
 586       implicit none
 587       character*(*) rekord,lancuch
 588       double precision wartosc,default
 589       integer ilen,iread
 590       external ilen
 591       iread=index(rekord,lancuch)
 592       if (iread.eq.0) then
 593         wartosc=default
 594         return
 595       endif
 596       iread=iread+ilen(lancuch)+1
 597       read (rekord(iread:),*) wartosc
 598       return
 599       end
 600 c----------------------------------------------------------------------------
 601       subroutine multreada(rekord,lancuch,tablica,dim,default)
 602       implicit none
 603       integer dim,i
 604       double precision tablica(dim),default
 605       character*(*) rekord,lancuch
 606       integer ilen,iread
 607       external ilen
 608       do i=1,dim
 609         tablica(i)=default
 610       enddo
 611       iread=index(rekord,lancuch)
 612       if (iread.eq.0) return
 613       iread=iread+ilen(lancuch)+1
 614       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 615    10 return
 616       end
 617 c----------------------------------------------------------------------------
 618       subroutine readi(rekord,lancuch,wartosc,default)
 619       implicit none
 620       character*(*) rekord,lancuch
 621       integer wartosc,default
 622       integer ilen,iread
 623       external ilen
 624       iread=index(rekord,lancuch)
 625       if (iread.eq.0) then
 626         wartosc=default
 627         return
 628       endif
 629       iread=iread+ilen(lancuch)+1
 630       read (rekord(iread:),*) wartosc
 631       return
 632       end
 633 c----------------------------------------------------------------------------
 634       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 635       implicit none
 636       integer dim,i
 637       integer tablica(dim),default
 638       character*(*) rekord,lancuch
 639       character*80 aux
 640       integer ilen,iread
 641       external ilen
 642       do i=1,dim
 643         tablica(i)=default
 644       enddo
 645       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 646       if (iread.eq.0) return
 647       iread=iread+ilen(lancuch)+1
 648       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 649    10 return
 650       end
 651 c----------------------------------------------------------------------------
 652       subroutine card_concat(card)
 653       include 'DIMENSIONS'
 654       include 'COMMON.IOUNITS'
 655       character*(*) card
 656       character*80 karta,ucase
 657       external ilen
 658       read (inp,'(a)') karta
 659       karta=ucase(karta)
 660       card=' '
 661       do while (karta(80:80).eq.'&')
 662         card=card(:ilen(card)+1)//karta(:79)
 663         read (inp,'(a)') karta
 664         karta=ucase(karta)
 665       enddo
 666       card=card(:ilen(card)+1)//karta
 667       return
 668       end
 669 c----------------------------------------------------------------------------
 670       subroutine openunits
 671       implicit none
 672       include 'DIMENSIONS'
 673 #ifdef MPI
 674       include "mpif.h"
 675       character*3 liczba
 676       include "COMMON.MPI"
 677 #endif
 678       include 'COMMON.IOUNITS'
 679       include 'COMMON.CONTROL'
 680       integer lenpre,lenpot,ilen
 681       external ilen
 682       character*16 cformat,cprint
 683       character*16 ucase
 684       integer lenint,lenout
 685       call getenv('INPUT',prefix)
 686       call getenv('OUTPUT',prefout)
 687       call getenv('INTIN',prefintin)
 688       call getenv('COORD',cformat)
 689       call getenv('PRINTCOOR',cprint)
 690       call getenv('SCRATCHDIR',scratchdir)
 691       from_bx=.true.
 692       from_cx=.false.
 693       if (index(ucase(cformat),'CX').gt.0) then
 694         from_cx=.true.
 695         from_bx=.false.
 696       endif
 697       from_cart=.true.
 698       lenpre=ilen(prefix)
 699       lenout=ilen(prefout)
 700       lenint=ilen(prefintin)
 701 C Get the names and open the input files
 702       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 703 #ifdef MPI
 704       write (liczba,'(bz,i3.3)') me
 705       outname=prefout(:lenout)//'_clust.out_'//liczba
 706 #else
 707       outname=prefout(:lenout)//'_clust.out'
 708 #endif
 709       if (from_bx) then
 710         intinname=prefintin(:lenint)//'.bx'
 711       else if (from_cx) then
 712         intinname=prefintin(:lenint)//'.cx'
 713       else
 714         intinname=prefintin(:lenint)//'.int'
 715       endif
 716       rmsname=prefintin(:lenint)//'.rms'
 717       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 718      &       status='unknown')
 719       open (jrms,file=rmsname,status='unknown')
 720       open(iout,file=outname,status='unknown')
 721 C Get parameter filenames and open the parameter files.
 722       call getenv('BONDPAR',bondname)
 723       open (ibond,file=bondname,status='old')
 724       call getenv('THETPAR',thetname)
 725       open (ithep,file=thetname,status='old')
 726       call getenv('ROTPAR',rotname)
 727       open (irotam,file=rotname,status='old')
 728       call getenv('TORPAR',torname)
 729       open (itorp,file=torname,status='old')
 730       call getenv('TORDPAR',tordname)
 731       open (itordp,file=tordname,status='old')
 732       call getenv('FOURIER',fouriername)
 733       open (ifourier,file=fouriername,status='old')
 734       call getenv('ELEPAR',elename)
 735       open (ielep,file=elename,status='old')
 736       call getenv('SIDEPAR',sidename)
 737       open (isidep,file=sidename,status='old')
 738       call getenv('SIDEP',sidepname)
 739       open (isidep1,file=sidepname,status="old")
 740       call getenv('SCCORPAR',sccorname)
 741       open (isccor,file=sccorname,status="old")
 742 #ifndef OLDSCP
 743 C
 744 C 8/9/01 In the newest version SCp interaction constants are read from a file
 745 C Use -DOLDSCP to use hard-coded constants instead.
 746 C
 747       call getenv('SCPPAR',scpname)
 748       open (iscpp,file=scpname,status='old')
 749 #endif
 750       return
 751       end
 752 c-------------------------------------------------------------------------------
 753       subroutine read_dist_constr
 754       implicit real*8 (a-h,o-z)
 755       include 'DIMENSIONS'
 756       include 'COMMON.CONTROL'
 757       include 'COMMON.CHAIN'
 758       include 'COMMON.IOUNITS'
 759       include 'COMMON.SBRIDGE'
 760       integer ifrag_(2,100),ipair_(2,100)
 761       double precision wfrag_(100),wpair_(100)
 762       character*500 controlcard
 763 c      write (iout,*) "Calling read_dist_constr"
 764 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 765       call card_concat(controlcard)
 766 c      call flush(iout)
 767 c      write (iout,'(a)') controlcard
 768       call readi(controlcard,"NFRAG",nfrag_,0)
 769       call readi(controlcard,"NPAIR",npair_,0)
 770       call readi(controlcard,"NDIST",ndist_,0)
 771       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 772       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 773       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 774       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 775       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 776       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 777       write (iout,*) "IFRAG"
 778       do i=1,nfrag_
 779         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 780       enddo
 781       write (iout,*) "IPAIR"
 782       do i=1,npair_
 783         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 784       enddo
 785       call flush(iout)
 786       if (.not.refstr .and. nfrag_.gt.0) then
 787         write (iout,*)
 788      &  "ERROR: no reference structure to compute distance restraints"
 789         write (iout,*)
 790      &  "Restraints must be specified explicitly (NDIST=number)"
 791         stop
 792       endif
 793       if (nfrag_.lt.2 .and. npair_.gt.0) then
 794         write (iout,*) "ERROR: Less than 2 fragments specified",
 795      &   " but distance restraints between pairs requested"
 796         stop
 797       endif
 798       call flush(iout)
 799       do i=1,nfrag_
 800         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 801         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 802      &    ifrag_(2,i)=nstart_sup+nsup-1
 803 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 804         call flush(iout)
 805         if (wfrag_(i).gt.0.0d0) then
 806         do j=ifrag_(1,i),ifrag_(2,i)-1
 807           do k=j+1,ifrag_(2,i)
 808             write (iout,*) "j",j," k",k
 809             ddjk=dist(j,k)
 810             if (constr_dist.eq.1) then
 811             nhpb=nhpb+1
 812             ihpb(nhpb)=j
 813             jhpb(nhpb)=k
 814               dhpb(nhpb)=ddjk
 815             forcon(nhpb)=wfrag_(i)
 816             else if (constr_dist.eq.2) then
 817               if (ddjk.le.dist_cut) then
 818                 nhpb=nhpb+1
 819                 ihpb(nhpb)=j
 820                 jhpb(nhpb)=k
 821                 dhpb(nhpb)=ddjk
 822                 forcon(nhpb)=wfrag_(i)
 823               endif
 824             else
 825               nhpb=nhpb+1
 826               ihpb(nhpb)=j
 827               jhpb(nhpb)=k
 828               dhpb(nhpb)=ddjk
 829               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 830             endif
 831             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 832      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 833           enddo
 834         enddo
 835         endif
 836       enddo
 837       do i=1,npair_
 838         if (wpair_(i).gt.0.0d0) then
 839         ii = ipair_(1,i)
 840         jj = ipair_(2,i)
 841         if (ii.gt.jj) then
 842           itemp=ii
 843           ii=jj
 844           jj=itemp
 845         endif
 846         do j=ifrag_(1,ii),ifrag_(2,ii)
 847           do k=ifrag_(1,jj),ifrag_(2,jj)
 848             nhpb=nhpb+1
 849             ihpb(nhpb)=j
 850             jhpb(nhpb)=k
 851             forcon(nhpb)=wpair_(i)
 852             dhpb(nhpb)=dist(j,k)
 853             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 854      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 855           enddo
 856         enddo
 857         endif
 858       enddo
 859       do i=1,ndist_
 860         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 861      &     ibecarb(i),forcon(nhpb+1)
 862         if (forcon(nhpb+1).gt.0.0d0) then
 863           nhpb=nhpb+1
 864           if (ibecarb(i).gt.0) then
 865             ihpb(i)=ihpb(i)+nres
 866             jhpb(i)=jhpb(i)+nres
 867           endif
 868           if (dhpb(nhpb).eq.0.0d0)
 869      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 870         endif
 871       enddo
 872       do i=1,nhpb
 873           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 874      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 875       enddo
 876       call flush(iout)
 877       return
 878       end
 879 <<<<<<< HEAD
 880 =======
 881
 882 c====-------------------------------------------------------------------
 883       subroutine read_constr_homology(lprn)
 884
 885       include 'DIMENSIONS'
 886 #ifdef MPI
 887       include 'mpif.h'
 888 #endif
 889       include 'COMMON.SETUP'
 890       include 'COMMON.CONTROL'
 891       include 'COMMON.CHAIN'
 892       include 'COMMON.IOUNITS'
 893       include 'COMMON.GEO'
 894       include 'COMMON.INTERACT'
 895       include 'COMMON.HOMRESTR'
 896 c
 897 c For new homol impl
 898 c
 899       include 'COMMON.VAR'
 900 c     include 'include_unres/COMMON.VAR'
 901 c
 902
 903 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 904 c    &                 dist_cut
 905 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 906 c    &    sigma_odl_temp(maxres,maxres,max_template)
 907       character*2 kic2
 908       character*24 model_ki_dist, model_ki_angle
 909       character*500 controlcard
 910       integer ki, i, j, k, l
 911       logical lprn
 912       logical unres_pdb
 913 c
 914 c     FP - Nov. 2014 Temporary specifications for new vars
 915 c
 916       double precision rescore_tmp,x12,y12,z12
 917       double precision, dimension (max_template,maxres) :: rescore
 918       character*24 tpl_k_rescore
 919 c -----------------------------------------------------------------
 920 c Reading multiple PDB ref structures and calculation of retraints
 921 c not using pre-computed ones stored in files model_ki_{dist,angle}
 922 c FP (Nov., 2014)
 923 c -----------------------------------------------------------------
 924 c
 925 c
 926 c Alternative: reading from input
 927 #ifdef DEBUG
 928       write (iout,*) "BEGIN READ HOMOLOGY INFO"
 929 #ifdef AIX
 930       call flush_(iout)
 931 #else
 932       call flush(iout)
 933 #endif
 934 #endif
 935       call card_concat(controlcard)
 936       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 937       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 938       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 939       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 940       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 941
 942       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
 943       if (homol_nset.gt.1)then
 944          call readi(controlcard,"ISET",iset,1)
 945          call card_concat(controlcard)
 946          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 947       else
 948         iset=1
 949         waga_homology(1)=1.0
 950       endif
 951 c
 952 #ifdef DEBUG
 953       write(iout,*) "read_constr_homology iset",iset
 954       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 955 #ifdef AIX
 956       call flush_(iout)
 957 #else
 958       call flush(iout)
 959 #endif
 960 #endif
 961 cd      write (iout,*) "nnt",nnt," nct",nct
 962 cd      call flush(iout)
 963
 964
 965       lim_odl=0
 966       lim_dih=0
 967 c
 968 c  New
 969 c
 970       lim_theta=0
 971       lim_xx=0
 972 c
 973 c  Reading HM global scores (prob not required)
 974 c
 975 c     open (4,file="HMscore")
 976 c     do k=1,constr_homology
 977 c       read (4,*,end=521) hmscore_tmp
 978 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 979 c       write(*,*) "Model", k, ":", hmscore(k)
 980 c     enddo
 981 c521  continue
 982
 983 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 984
 985       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
 986       do k=1,constr_homology
 987
 988         read(inp,'(a)') pdbfile
 989 c        write (iout,*) "k ",k," pdbfile ",pdbfile
 990 c  Next stament causes error upon compilation (?)
 991 c       if(me.eq.king.or. .not. out1file)
 992 c         write (iout,'(2a)') 'PDB data will be read from file ',
 993 c    &   pdbfile(:ilen(pdbfile))
 994         open(ipdbin,file=pdbfile,status='old',err=33)
 995         goto 34
 996   33    write (iout,'(a)') 'Error opening PDB file.'
 997         stop
 998   34    continue
 999 c        print *,'Begin reading pdb data'
1000 c
1001 c Files containing res sim or local scores (former containing sigmas)
1002 c
1003
1004         write(kic2,'(bz,i2.2)') k
1005
1006         tpl_k_rescore="template"//kic2//".sco"
1007 c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
1008 c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
1009 c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
1010 c       tpl_k_sigma_d="template"//kic2//".sigma_d"
1011
1012         unres_pdb=.false.
1013         call readpdb
1014         do i=1,2*nres
1015           do j=1,3
1016             crefjlee(j,i)=c(j,i)
1017           enddo
1018         enddo
1019 #ifdef DEBUG
1020         do i=1,nres
1021           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1022      &      (crefjlee(j,i+nres),j=1,3)
1023         enddo
1024         write (iout,*) "READ HOMOLOGY INFO"
1025         write (iout,*) "read_constr_homology x: after reading pdb file"
1026         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1027         write (iout,*) "waga_dist",waga_dist
1028         write (iout,*) "waga_angle",waga_angle
1029         write (iout,*) "waga_theta",waga_theta
1030         write (iout,*) "waga_d",waga_d
1031         write (iout,*) "dist_cut",dist_cut
1032 #endif
1033 #ifdef AIX
1034       call flush_(iout)
1035 #else
1036       call flush(iout)
1037 #endif
1038
1039 c
1040 c     Distance restraints
1041 c
1042 c          ... --> odl(k,ii)
1043 C Copy the coordinates from reference coordinates (?)
1044         do i=1,2*nres
1045           do j=1,3
1046             c(j,i)=cref(j,i)
1047 c           write (iout,*) "c(",j,i,") =",c(j,i)
1048           enddo
1049         enddo
1050 c
1051 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1052 c
1053 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1054           open (ientin,file=tpl_k_rescore,status='old')
1055           do irec=1,maxdim ! loop for reading res sim
1056             if (irec.eq.1) then
1057                rescore(k,irec)=0.0d0
1058                goto 1301
1059             endif
1060             read (ientin,*,end=1401) rescore_tmp
1061 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1062             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1063 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1064  1301     continue
1065           enddo
1066  1401   continue
1067           close (ientin)
1068 c         open (ientin,file=tpl_k_sigma_odl,status='old')
1069 c         do irec=1,maxdim ! loop for reading sigma_odl
1070 c            read (ientin,*,end=1401) i, j,
1071 c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
1072 c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
1073 c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k)
1074 c         enddo
1075 c 1401   continue
1076 c         close (ientin)
1077         if (waga_dist.ne.0.0d0) then
1078           ii=0
1079           do i = nnt,nct-2 ! right? without parallel.
1080             do j=i+2,nct ! right?
1081 c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
1082 c           do j=i+2,nres ! ibid
1083 c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
1084 c           do j=i+2,nct ! ibid
1085               ii=ii+1
1086 c             write (iout,*) "k",k
1087 c             write (iout,*) "i",i," j",j," constr_homology",
1088 c    &                       constr_homology
1089               ires_homo(ii)=i
1090               jres_homo(ii)=j
1091 c
1092 c Attempt to replace dist(i,j) by its definition in ...
1093 c
1094               x12=c(1,i)-c(1,j)
1095               y12=c(2,i)-c(2,j)
1096               z12=c(3,i)-c(3,j)
1097               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1098               odl(k,ii)=distal
1099 c
1100 c             odl(k,ii)=dist(i,j)
1101 c             write (iout,*) "dist(",i,j,") =",dist(i,j)
1102 c             write (iout,*) "distal = ",distal
1103 c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
1104 c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1105 c    &                      "rescore(",k,j,") =",rescore(k,j)
1106 c
1107 c  Calculation of sigma from res sim
1108 c
1109 c             if (odl(k,ii).le.6.0d0) then
1110 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
1111 c    Other functional forms possible depending on odl(k,ii), eg.
1112 c
1113             if (odl(k,ii).le.dist_cut) then
1114               sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
1115 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
1116             else
1117 #ifdef OLDSIGMA
1118               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
1119      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1120 #else
1121               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
1122      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1123 #endif
1124 c   Following expr replaced by a positive exp argument
1125 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1126 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
1127
1128 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
1129 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
1130             endif
1131 c
1132               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
1133 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
1134 c
1135 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
1136 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
1137             enddo
1138 c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
1139           enddo
1140 c         lim_odl=ii
1141 c         if (constr_homology.gt.0) call homology_partition
1142         endif
1143 c
1144 c     Theta, dihedral and SC retraints
1145 c
1146         if (waga_angle.gt.0.0d0) then
1147 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1148 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1149 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1150 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1151 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1152 c    &                            sigma_dih(k,i+nnt-1)
1153 c         enddo
1154 c1402   continue
1155 c         close (ientin)
1156           do i = nnt+3,nct ! right? without parallel.
1157 c         do i=1,nres ! alternative for bounds acc to readpdb?
1158 c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
1159 c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
1160             dih(k,i)=phiref(i) ! right?
1161 c           read (ientin,*) sigma_dih(k,i) ! original variant
1162 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1163 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1164 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1165 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1166 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1167
1168             sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
1169      &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
1170 c
1171 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1172 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1173 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1174 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1175 c   Instead of res sim other local measure of b/b str reliability possible
1176             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1177 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1178             if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
1179 c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
1180           enddo
1181         endif
1182
1183         if (waga_theta.gt.0.0d0) then
1184 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1185 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1186 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1187 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1188 c    &                              sigma_theta(k,i+nnt-1)
1189 c         enddo
1190 c1403   continue
1191 c         close (ientin)
1192
1193           do i = nnt+2,nct ! right? without parallel.
1194 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1195 c         do i=ithet_start,ithet_end ! with FG parallel.
1196              thetatpl(k,i)=thetaref(i)
1197 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1198 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1199 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1200 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1201 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1202              sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
1203      &                        rescore(k,i-2) !  right expression ?
1204              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1205
1206 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1207 c                             rescore(k,i-2) !  right expression ?
1208 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1209              if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
1210           enddo
1211         endif
1212
1213         if (waga_d.gt.0.0d0) then
1214 c       open (ientin,file=tpl_k_sigma_d,status='old')
1215 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1216 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1217 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1218 c    &                          sigma_d(k,i+nnt-1)
1219 c         enddo
1220 c1404   continue
1221           close (ientin)
1222
1223           do i = nnt,nct ! right? without parallel.
1224 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1225 c         do i=loc_start,loc_end ! with FG parallel.
1226              if (itype(i).eq.10) goto 1 ! right?
1227                xxtpl(k,i)=xxref(i)
1228                yytpl(k,i)=yyref(i)
1229                zztpl(k,i)=zzref(i)
1230 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1231 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1232 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1233 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1234                sigma_d(k,i)=rescore(k,i) !  right expression ?
1235                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1236
1237 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1238 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1239 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1240                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1241     1     continue
1242           enddo
1243         endif
1244         close(ientin)
1245       enddo
1246       if (waga_dist.ne.0.0d0) lim_odl=ii
1247       if (constr_homology.gt.0) call homology_partition
1248       if (constr_homology.gt.0) call init_int_table
1249 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1250 cd     & "lim_xx=",lim_xx
1251 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1252 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1253 c
1254 c Print restraints
1255 c
1256       if (.not.lprn) return
1257 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1258       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1259        write (iout,*) "Distance restraints from templates"
1260        do ii=1,lim_odl
1261        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
1262      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
1263        enddo
1264        write (iout,*) "Dihedral angle restraints from templates"
1265        do i=nnt+3,lim_dih
1266         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
1267      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1268        enddo
1269        write (iout,*) "Virtual-bond angle restraints from templates"
1270        do i=nnt+2,lim_theta
1271         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
1272      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1273        enddo
1274        write (iout,*) "SC restraints from templates"
1275        do i=nnt,lim_xx
1276         write(iout,'(i5,10(4f8.2,4x))') i,
1277      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1278      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1279        enddo
1280       endif
1281 c -----------------------------------------------------------------
1282       return
1283       end
1284
1285
1286 >>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit