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