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