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