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