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