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       logical normalize,next
 851       integer restr_type
 852       double precision xlink(4,0:4) /
 853 c           a          b       c     sigma
 854      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
 855      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
 856      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
 857      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
 858      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
 859       write (iout,*) "Calling read_dist_constr"
 860 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 861 c      call flush(iout)
 862       next=.true.
 863
 864       DO WHILE (next)
 865
 866       call card_concat(controlcard)
 867       next = index(controlcard,"NEXT").gt.0
 868       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
 869       write (iout,*) "restr_type",restr_type
 870       call readi(controlcard,"NFRAG",nfrag_,0)
 871       call readi(controlcard,"NFRAG",nfrag_,0)
 872       call readi(controlcard,"NPAIR",npair_,0)
 873       call readi(controlcard,"NDIST",ndist_,0)
 874       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 875       if (restr_type.eq.10)
 876      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
 877       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 878       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 879       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 880       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 881       normalize = index(controlcard,"NORMALIZE").gt.0
 882       write (iout,*) "WBOLTZD",wboltzd
 883 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 884 c      write (iout,*) "IFRAG"
 885 c      do i=1,nfrag_
 886 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 887 c      enddo
 888 c      write (iout,*) "IPAIR"
 889 c      do i=1,npair_
 890 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 891 c      enddo
 892       if (nfrag_.gt.0) then
 893         nres0=nres
 894         read(inp,'(a)') pdbfile
 895         write (iout,*)
 896      & "Distance restraints will be constructed from structure ",pdbfile
 897         open(ipdbin,file=pdbfile,status='old',err=11)
 898         call readpdb(.true.)
 899         nres=nres0
 900         close(ipdbin)
 901       endif
 902       do i=1,nfrag_
 903         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 904         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 905      &    ifrag_(2,i)=nstart_sup+nsup-1
 906 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 907 c        call flush(iout)
 908         if (wfrag_(i).eq.0.0d0) cycle
 909         do j=ifrag_(1,i),ifrag_(2,i)-1
 910           do k=j+1,ifrag_(2,i)
 911 c            write (iout,*) "j",j," k",k
 912             ddjk=dist(j,k)
 913             if (restr_type.eq.1) then
 914               nhpb=nhpb+1
 915               irestr_type(nhpb)=1
 916               ihpb(nhpb)=j
 917               jhpb(nhpb)=k
 918               dhpb(nhpb)=ddjk
 919               forcon(nhpb)=wfrag_(i)
 920             else if (constr_dist.eq.2) then
 921               if (ddjk.le.dist_cut) then
 922                 nhpb=nhpb+1
 923                 irestr_type(nhpb)=1
 924                 ihpb(nhpb)=j
 925                 jhpb(nhpb)=k
 926                 dhpb(nhpb)=ddjk
 927                 forcon(nhpb)=wfrag_(i)
 928               endif
 929             else if (restr_type.eq.3) then
 930               nhpb=nhpb+1
 931               irestr_type(nhpb)=1
 932               ihpb(nhpb)=j
 933               jhpb(nhpb)=k
 934               dhpb(nhpb)=ddjk
 935               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 936             endif
 937             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
 938      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 939           enddo
 940         enddo
 941       enddo
 942       do i=1,npair_
 943         if (wpair_(i).eq.0.0d0) cycle
 944         ii = ipair_(1,i)
 945         jj = ipair_(2,i)
 946         if (ii.gt.jj) then
 947           itemp=ii
 948           ii=jj
 949           jj=itemp
 950         endif
 951         do j=ifrag_(1,ii),ifrag_(2,ii)
 952           do k=ifrag_(1,jj),ifrag_(2,jj)
 953             if (restr_type.eq.1) then
 954               nhpb=nhpb+1
 955               irestr_type(nhpb)=1
 956               ihpb(nhpb)=j
 957               jhpb(nhpb)=k
 958               dhpb(nhpb)=ddjk
 959               forcon(nhpb)=wfrag_(i)
 960             else if (constr_dist.eq.2) then
 961               if (ddjk.le.dist_cut) then
 962                 nhpb=nhpb+1
 963                 irestr_type(nhpb)=1
 964                 ihpb(nhpb)=j
 965                 jhpb(nhpb)=k
 966                 dhpb(nhpb)=ddjk
 967                 forcon(nhpb)=wfrag_(i)
 968               endif
 969             else if (restr_type.eq.3) then
 970               nhpb=nhpb+1
 971               irestr_type(nhpb)=1
 972               ihpb(nhpb)=j
 973               jhpb(nhpb)=k
 974               dhpb(nhpb)=ddjk
 975               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 976             endif
 977             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
 978      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 979           enddo
 980         enddo
 981       enddo
 982
 983 c      print *,ndist_
 984       write (iout,*) "Distance restraints as read from input"
 985       do i=1,ndist_
 986         if (restr_type.eq.11) then
 987           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
 988      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 989 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 990           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
 991           nhpb=nhpb+1
 992           irestr_type(nhpb)=11
 993           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
 994      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 995      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 996           if (ibecarb(nhpb).gt.0) then
 997             ihpb(nhpb)=ihpb(nhpb)+nres
 998             jhpb(nhpb)=jhpb(nhpb)+nres
 999           endif
1000         else if (constr_dist.eq.10) then
1001 c Cross-lonk Markov-like potential
1002           call card_concat(controlcard)
1003           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1004           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1005           ibecarb(nhpb+1)=0
1006           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1007           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1008           if (index(controlcard,"ZL").gt.0) then
1009             link_type=1
1010           else if (index(controlcard,"ADH").gt.0) then
1011             link_type=2
1012           else if (index(controlcard,"PDH").gt.0) then
1013             link_type=3
1014           else if (index(controlcard,"DSS").gt.0) then
1015             link_type=4
1016           else
1017             link_type=0
1018           endif
1019           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1020      &       xlink(1,link_type))
1021           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1022      &       xlink(2,link_type))
1023           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1024      &       xlink(3,link_type))
1025           call reada(controlcard,"SIGMA",forcon(nhpb+1),
1026      &       xlink(4,link_type))
1027           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1028 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1029 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1030           if (forcon(nhpb+1).le.0.0d0 .or.
1031      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1032           nhpb=nhpb+1
1033           irestr_type(nhpb)=10
1034           if (ibecarb(nhpb).gt.0) then
1035             ihpb(nhpb)=ihpb(nhpb)+nres
1036             jhpb(nhpb)=jhpb(nhpb)+nres
1037           endif
1038           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1039      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1040      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1041      &     irestr_type(nhpb)
1042         else
1043 C        print *,"in else"
1044           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1045      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1046           if (forcon(nhpb+1).gt.0.0d0) then
1047           nhpb=nhpb+1
1048           if (dhpb1(nhpb).eq.0.0d0) then
1049             irestr_type(nhpb)=1
1050           else
1051             irestr_type(nhpb)=2
1052           endif
1053           if (ibecarb(nhpb).gt.0) then
1054             ihpb(nhpb)=ihpb(nhpb)+nres
1055             jhpb(nhpb)=jhpb(nhpb)+nres
1056           endif
1057           if (dhpb(nhpb).eq.0.0d0)
1058      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1059         endif
1060         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1061      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1062         endif
1063 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1064 C        if (forcon(nhpb+1).gt.0.0d0) then
1065 C          nhpb=nhpb+1
1066 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1067       enddo
1068
1069       ENDDO ! next
1070
1071       fordepthmax=0.0d0
1072       if (normalize) then
1073         do i=nss+1,nhpb
1074           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1075      &      fordepthmax=fordepth(i)
1076         enddo
1077         do i=nss+1,nhpb
1078           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1079         enddo
1080       endif
1081       if (nhpb.gt.nss)  then
1082         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1083      &  "The following",nhpb-nss,
1084      &  " distance restraints have been imposed:",
1085      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
1086      &  "  score"," type"
1087         do i=nss+1,nhpb
1088           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1089      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1090      &  irestr_type(i)
1091         enddo
1092       endif
1093       call hpb_partition
1094       call flush(iout)
1095       return
1096    11 write (iout,*)"read_dist_restr: error reading reference structure"
1097       stop
1098       end
1099 c-------------------------------------------------------------------------------
1100       subroutine read_saxs_constr
1101       implicit real*8 (a-h,o-z)
1102       include 'DIMENSIONS'
1103 #ifdef MPI
1104       include 'mpif.h'
1105 #endif
1106       include 'COMMON.CONTROL'
1107       include 'COMMON.CHAIN'
1108       include 'COMMON.IOUNITS'
1109       include 'COMMON.SBRIDGE'
1110       double precision cm(3)
1111 c      read(inp,*) nsaxs
1112       write (iout,*) "Calling read_saxs nsaxs",nsaxs
1113       call flush(iout)
1114       if (saxs_mode.eq.0) then
1115 c SAXS distance distribution
1116       do i=1,nsaxs
1117         read(inp,*) distsaxs(i),Psaxs(i)
1118       enddo
1119       Cnorm = 0.0d0
1120       do i=1,nsaxs
1121         Cnorm = Cnorm + Psaxs(i)
1122       enddo
1123       write (iout,*) "Cnorm",Cnorm
1124       do i=1,nsaxs
1125         Psaxs(i)=Psaxs(i)/Cnorm
1126       enddo
1127       write (iout,*) "Normalized distance distribution from SAXS"
1128       do i=1,nsaxs
1129         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1130       enddo
1131       Wsaxs0=0.0d0
1132       do i=1,nsaxs
1133         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1134       enddo
1135       write (iout,*) "Wsaxs0",Wsaxs0
1136       else
1137 c SAXS "spheres".
1138       do i=1,nsaxs
1139         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1140       enddo
1141       do j=1,3
1142         cm(j)=0.0d0
1143       enddo
1144       do i=1,nsaxs
1145         do j=1,3
1146           cm(j)=cm(j)+Csaxs(j,i)
1147         enddo
1148       enddo
1149       do j=1,3
1150         cm(j)=cm(j)/nsaxs
1151       enddo
1152       do i=1,nsaxs
1153         do j=1,3
1154           Csaxs(j,i)=Csaxs(j,i)-cm(j)
1155         enddo
1156       enddo
1157       write (iout,*) "SAXS sphere coordinates"
1158       do i=1,nsaxs
1159         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
1160       enddo
1161       endif
1162       return
1163       end
1164 c====-------------------------------------------------------------------
1165       subroutine read_constr_homology(lprn)
1166
1167       include 'DIMENSIONS'
1168 #ifdef MPI
1169       include 'mpif.h'
1170 #endif
1171       include 'COMMON.SETUP'
1172       include 'COMMON.CONTROL'
1173       include 'COMMON.CHAIN'
1174       include 'COMMON.IOUNITS'
1175       include 'COMMON.GEO'
1176       include 'COMMON.INTERACT'
1177       include 'COMMON.NAMES'
1178       include 'COMMON.HOMRESTR'
1179 c
1180 c For new homol impl
1181 c
1182       include 'COMMON.VAR'
1183 c     include 'include_unres/COMMON.VAR'
1184 c
1185
1186 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1187 c    &                 dist_cut
1188 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1189 c    &    sigma_odl_temp(maxres,maxres,max_template)
1190       character*2 kic2
1191       character*24 model_ki_dist, model_ki_angle
1192       character*500 controlcard
1193       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1194       integer idomain(max_template,maxres)
1195       integer ilen
1196       external ilen
1197       logical lprn
1198       logical unres_pdb,liiflag
1199 c
1200 c     FP - Nov. 2014 Temporary specifications for new vars
1201 c
1202       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
1203      &                 rescore3_tmp
1204       double precision, dimension (max_template,maxres) :: rescore
1205       double precision, dimension (max_template,maxres) :: rescore2
1206       double precision, dimension (max_template,maxres) :: rescore3
1207       character*24 tpl_k_rescore
1208 c -----------------------------------------------------------------
1209 c Reading multiple PDB ref structures and calculation of retraints
1210 c not using pre-computed ones stored in files model_ki_{dist,angle}
1211 c FP (Nov., 2014)
1212 c -----------------------------------------------------------------
1213 c
1214 c
1215 c Alternative: reading from input
1216 #ifdef DEBUG
1217       write (iout,*) "BEGIN READ HOMOLOGY INFO"
1218 #ifdef AIX
1219       call flush_(iout)
1220 #else
1221       call flush(iout)
1222 #endif
1223 #endif
1224       call card_concat(controlcard)
1225       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1226       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1227       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1228       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1229       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1230       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1231       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1232       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1233       if (homol_nset.gt.1)then
1234          call readi(controlcard,"ISET",iset,1)
1235          call card_concat(controlcard)
1236          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1237       else
1238         iset=1
1239         waga_homology(1)=1.0
1240       endif
1241 c
1242 #ifdef DEBUG
1243       write(iout,*) "read_constr_homology iset",iset
1244       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1245 #ifdef AIX
1246       call flush_(iout)
1247 #else
1248       call flush(iout)
1249 #endif
1250 #endif
1251 cd      write (iout,*) "nnt",nnt," nct",nct
1252 cd      call flush(iout)
1253
1254
1255       lim_odl=0
1256       lim_dih=0
1257 c
1258 c  New
1259 c
1260 c
1261 c  Reading HM global scores (prob not required)
1262 c
1263       do i = nnt,nct
1264         do k=1,constr_homology
1265          idomain(k,i)=0
1266         enddo
1267       enddo
1268 c     open (4,file="HMscore")
1269 c     do k=1,constr_homology
1270 c       read (4,*,end=521) hmscore_tmp
1271 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1272 c       write(*,*) "Model", k, ":", hmscore(k)
1273 c     enddo
1274 c521  continue
1275
1276       ii=0
1277       do i = nnt,nct-2
1278         do j=i+2,nct
1279         ii=ii+1
1280         ii_in_use(ii)=0
1281         enddo
1282       enddo
1283 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1284
1285       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1286       do k=1,constr_homology
1287
1288         read(inp,'(a)') pdbfile
1289 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1290 c  Next stament causes error upon compilation (?)
1291 c       if(me.eq.king.or. .not. out1file)
1292 c         write (iout,'(2a)') 'PDB data will be read from file ',
1293 c    &   pdbfile(:ilen(pdbfile))
1294          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1295      &  pdbfile(:ilen(pdbfile))
1296         open(ipdbin,file=pdbfile,status='old',err=33)
1297         goto 34
1298   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1299      &  pdbfile(:ilen(pdbfile))
1300         stop
1301   34    continue
1302 c        print *,'Begin reading pdb data'
1303 c
1304 c Files containing res sim or local scores (former containing sigmas)
1305 c
1306
1307         write(kic2,'(bz,i2.2)') k
1308
1309         tpl_k_rescore="template"//kic2//".sco"
1310
1311         unres_pdb=.false.
1312         call readpdb_template(k)
1313         do i=1,2*nres
1314           do j=1,3
1315             crefjlee(j,i)=c(j,i)
1316           enddo
1317         enddo
1318 #ifdef DEBUG
1319         do i=1,nres
1320           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1321      &      (crefjlee(j,i+nres),j=1,3)
1322         enddo
1323         write (iout,*) "READ HOMOLOGY INFO"
1324         write (iout,*) "read_constr_homology x: after reading pdb file"
1325         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1326         write (iout,*) "waga_dist",waga_dist
1327         write (iout,*) "waga_angle",waga_angle
1328         write (iout,*) "waga_theta",waga_theta
1329         write (iout,*) "waga_d",waga_d
1330         write (iout,*) "dist_cut",dist_cut
1331 #endif
1332 #ifdef AIX
1333       call flush_(iout)
1334 #else
1335       call flush(iout)
1336 #endif
1337
1338 c
1339 c     Distance restraints
1340 c
1341 c          ... --> odl(k,ii)
1342 C Copy the coordinates from reference coordinates (?)
1343         do i=1,2*nres
1344           do j=1,3
1345 c            c(j,i)=cref(j,i)
1346 c           write (iout,*) "c(",j,i,") =",c(j,i)
1347           enddo
1348         enddo
1349 c
1350 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1351 c
1352 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1353           open (ientin,file=tpl_k_rescore,status='old')
1354           if (nnt.gt.1) rescore(k,1)=0.0d0
1355           do irec=nnt,nct ! loop for reading res sim
1356             if (read2sigma) then
1357              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1358      &                                rescore3_tmp,idomain_tmp
1359              i_tmp=i_tmp+nnt-1
1360              idomain(k,i_tmp)=idomain_tmp
1361              rescore(k,i_tmp)=rescore_tmp
1362              rescore2(k,i_tmp)=rescore2_tmp
1363              rescore3(k,i_tmp)=rescore3_tmp
1364              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1365      &                      i_tmp,rescore2_tmp,rescore_tmp,
1366      &                                rescore3_tmp,idomain_tmp
1367             else
1368              idomain(k,irec)=1
1369              read (ientin,*,end=1401) rescore_tmp
1370
1371 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1372              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1373 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1374             endif
1375           enddo
1376  1401   continue
1377         close (ientin)
1378         if (waga_dist.ne.0.0d0) then
1379           ii=0
1380           do i = nnt,nct-2
1381             do j=i+2,nct
1382
1383               x12=c(1,i)-c(1,j)
1384               y12=c(2,i)-c(2,j)
1385               z12=c(3,i)-c(3,j)
1386               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1387 c              write (iout,*) k,i,j,distal,dist2_cut
1388
1389             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1390      &            .and. distal.le.dist2_cut ) then
1391
1392               ii=ii+1
1393               ii_in_use(ii)=1
1394               l_homo(k,ii)=.true.
1395
1396 c             write (iout,*) "k",k
1397 c             write (iout,*) "i",i," j",j," constr_homology",
1398 c    &                       constr_homology
1399               ires_homo(ii)=i
1400               jres_homo(ii)=j
1401               odl(k,ii)=distal
1402               if (read2sigma) then
1403                 sigma_odl(k,ii)=0
1404                 do ik=i,j
1405                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1406                 enddo
1407                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1408                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1409      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1410               else
1411                 if (odl(k,ii).le.dist_cut) then
1412                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1413                 else
1414 #ifdef OLDSIGMA
1415                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1416      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1417 #else
1418                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1419      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1420 #endif
1421                 endif
1422               endif
1423               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1424             else
1425               ii=ii+1
1426               l_homo(k,ii)=.false.
1427             endif
1428             enddo
1429           enddo
1430         lim_odl=ii
1431         endif
1432 c
1433 c     Theta, dihedral and SC retraints
1434 c
1435         if (waga_angle.gt.0.0d0) then
1436 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1437 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1438 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1439 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1440 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1441 c    &                            sigma_dih(k,i+nnt-1)
1442 c         enddo
1443 c1402   continue
1444 c         close (ientin)
1445           do i = nnt+3,nct
1446             if (idomain(k,i).eq.0) then
1447                sigma_dih(k,i)=0.0
1448                cycle
1449             endif
1450             dih(k,i)=phiref(i) ! right?
1451 c           read (ientin,*) sigma_dih(k,i) ! original variant
1452 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1453 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1454 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1455 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1456 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1457
1458             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1459      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1460 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1461 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1462 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1463 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1464 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1465 c   Instead of res sim other local measure of b/b str reliability possible
1466             if (sigma_dih(k,i).ne.0)
1467      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1468 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1469           enddo
1470           lim_dih=nct-nnt-2
1471         endif
1472
1473         if (waga_theta.gt.0.0d0) then
1474 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1475 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1476 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1477 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1478 c    &                              sigma_theta(k,i+nnt-1)
1479 c         enddo
1480 c1403   continue
1481 c         close (ientin)
1482
1483           do i = nnt+2,nct ! right? without parallel.
1484 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1485 c         do i=ithet_start,ithet_end ! with FG parallel.
1486              if (idomain(k,i).eq.0) then
1487               sigma_theta(k,i)=0.0
1488               cycle
1489              endif
1490              thetatpl(k,i)=thetaref(i)
1491 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1492 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1493 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1494 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1495 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1496              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1497      &                        rescore(k,i-2))/3.0
1498 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1499              if (sigma_theta(k,i).ne.0)
1500      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1501
1502 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1503 c                             rescore(k,i-2) !  right expression ?
1504 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1505           enddo
1506         endif
1507
1508         if (waga_d.gt.0.0d0) then
1509 c       open (ientin,file=tpl_k_sigma_d,status='old')
1510 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1511 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1512 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1513 c    &                          sigma_d(k,i+nnt-1)
1514 c         enddo
1515 c1404   continue
1516
1517           do i = nnt,nct ! right? without parallel.
1518 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1519 c         do i=loc_start,loc_end ! with FG parallel.
1520                if (itype(i).eq.10) cycle
1521                if (idomain(k,i).eq.0 ) then
1522                   sigma_d(k,i)=0.0
1523                   cycle
1524                endif
1525                xxtpl(k,i)=xxref(i)
1526                yytpl(k,i)=yyref(i)
1527                zztpl(k,i)=zzref(i)
1528 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1529 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1530 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1531 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1532                sigma_d(k,i)=rescore3(k,i) !  right expression ?
1533                if (sigma_d(k,i).ne.0)
1534      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1535
1536 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1537 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1538 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1539           enddo
1540         endif
1541       enddo
1542 c
1543 c remove distance restraints not used in any model from the list
1544 c shift data in all arrays
1545 c
1546       if (waga_dist.ne.0.0d0) then
1547         ii=0
1548         liiflag=.true.
1549         do i=nnt,nct-2
1550          do j=i+2,nct
1551           ii=ii+1
1552           if (ii_in_use(ii).eq.0.and.liiflag) then
1553             liiflag=.false.
1554             iistart=ii
1555           endif
1556           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1557      &                   .not.liiflag.and.ii.eq.lim_odl) then
1558              if (ii.eq.lim_odl) then
1559               iishift=ii-iistart+1
1560              else
1561               iishift=ii-iistart
1562              endif
1563              liiflag=.true.
1564              do ki=iistart,lim_odl-iishift
1565               ires_homo(ki)=ires_homo(ki+iishift)
1566               jres_homo(ki)=jres_homo(ki+iishift)
1567               ii_in_use(ki)=ii_in_use(ki+iishift)
1568               do k=1,constr_homology
1569                odl(k,ki)=odl(k,ki+iishift)
1570                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1571                l_homo(k,ki)=l_homo(k,ki+iishift)
1572               enddo
1573              enddo
1574              ii=ii-iishift
1575              lim_odl=lim_odl-iishift
1576           endif
1577          enddo
1578         enddo
1579       endif
1580       if (constr_homology.gt.0) call homology_partition
1581       if (constr_homology.gt.0) call init_int_table
1582 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1583 cd     & "lim_xx=",lim_xx
1584 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1585 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1586 c
1587 c Print restraints
1588 c
1589       if (.not.lprn) return
1590 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1591       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1592        write (iout,*) "Distance restraints from templates"
1593        do ii=1,lim_odl
1594        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1595      &  ii,ires_homo(ii),jres_homo(ii),
1596      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1597      &  ki=1,constr_homology)
1598        enddo
1599        write (iout,*) "Dihedral angle restraints from templates"
1600        do i=nnt+3,nct
1601         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1602      &      (rad2deg*dih(ki,i),
1603      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1604        enddo
1605        write (iout,*) "Virtual-bond angle restraints from templates"
1606        do i=nnt+2,nct
1607         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1608      &      (rad2deg*thetatpl(ki,i),
1609      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1610        enddo
1611        write (iout,*) "SC restraints from templates"
1612        do i=nnt,nct
1613         write(iout,'(i5,100(4f8.2,4x))') i,
1614      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1615      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1616        enddo
1617       endif
1618 c -----------------------------------------------------------------
1619       return
1620       end