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