source/cluster/wham/src-restraint/readrtns.F

   1       subroutine read_control
   2 C
   3 C Read molecular data
   4 C
   5       implicit none
   6       include 'DIMENSIONS'
   7       include 'sizesclu.dat'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.TIME1'
  10       include 'COMMON.SBRIDGE'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.HEADER'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.FREE'
  17       character*320 controlcard,ucase
  18 #ifdef MPL
  19       include 'COMMON.INFO'
  20 #endif
  21       integer i
  22
  23       read (INP,'(a80)') titel
  24       call card_concat(controlcard)
  25
  26       call readi(controlcard,'NRES',nres,0)
  27       call readi(controlcard,'RESCALE',rescale_mode,2)
  28       call readi(controlcard,'PDBOUT',outpdb,0)
  29       call readi(controlcard,'MOL2OUT',outmol2,0)
  30       refstr=(index(controlcard,'REFSTR').gt.0)
  31       write (iout,*) "REFSTR",refstr
  32       pdbref=(index(controlcard,'PDBREF').gt.0)
  33       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
  34       iscode=index(controlcard,'ONE_LETTER')
  35       tree=(index(controlcard,'MAKE_TREE').gt.0)
  36       min_var=(index(controlcard,'MINVAR').gt.0)
  37       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
  38       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
  39       call readi(controlcard,'NCUT',ncut,1)
  40       call readi(controlcard,'NSTART',nstart,0)
  41       call readi(controlcard,'NEND',nend,0)
  42       call reada(controlcard,'ECUT',ecut,10.0d0)
  43       call reada(controlcard,'PROB',prob_limit,0.99d0)
  44       write (iout,*) "Probability limit",prob_limit
  45       lgrp=(index(controlcard,'LGRP').gt.0)
  46       caonly=(index(controlcard,'CA_ONLY').gt.0)
  47       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
  48       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
  49       call readi(controlcard,'IOPT',iopt,2)
  50       lside = index(controlcard,"SIDE").gt.0
  51       efree = index(controlcard,"EFREE").gt.0
  52       call readi(controlcard,'NTEMP',nT,1)
  53       write (iout,*) "nT",nT
  54       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
  55       write (iout,*) "nT",nT
  56       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
  57       do i=1,nT
  58         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
  59       enddo
  60       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
  61       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
  62       lprint_int=index(controlcard,"PRINT_INT") .gt.0
  63       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  64       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  65       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  66      & " CONSTR_DIST",constr_dist
  67
  68       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
  69       write (iout,*) "with_homology_constr ",with_dihed_constr,
  70      & " CONSTR_HOMOLOGY",constr_homology
  71
  72       call flush(iout)
  73       if (min_var) iopt=1
  74       return
  75       end
  76 c--------------------------------------------------------------------------
  77       subroutine molread
  78 C
  79 C Read molecular data.
  80 C
  81       implicit none
  82       include 'DIMENSIONS'
  83       include 'COMMON.IOUNITS'
  84       include 'COMMON.GEO'
  85       include 'COMMON.VAR'
  86       include 'COMMON.INTERACT'
  87       include 'COMMON.LOCAL'
  88       include 'COMMON.NAMES'
  89       include 'COMMON.CHAIN'
  90       include 'COMMON.FFIELD'
  91       include 'COMMON.SBRIDGE'
  92       include 'COMMON.HEADER'
  93       include 'COMMON.CONTROL'
  94       include 'COMMON.CONTACTS'
  95       include 'COMMON.TIME1'
  96       include 'COMMON.TORCNSTR'
  97 #ifdef MPL
  98       include 'COMMON.INFO'
  99 #endif
 100       character*4 sequence(maxres)
 101       character*800 weightcard
 102       integer rescode
 103       double precision x(maxvar)
 104       integer itype_pdb(maxres)
 105       logical seq_comp
 106       integer i,j
 107 C
 108 C Body
 109 C
 110 C Read weights of the subsequent energy terms.
 111       call card_concat(weightcard)
 112       call reada(weightcard,'WSC',wsc,1.0d0)
 113       call reada(weightcard,'WLONG',wsc,wsc)
 114       call reada(weightcard,'WSCP',wscp,1.0d0)
 115       call reada(weightcard,'WELEC',welec,1.0D0)
 116       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 117       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 118       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 119       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 120       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 121       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 122       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 123       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 124       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 125       call reada(weightcard,'WBOND',wbond,1.0D0)
 126       call reada(weightcard,'WTOR',wtor,1.0D0)
 127       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 128       call reada(weightcard,'WANG',wang,1.0D0)
 129       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 130       call reada(weightcard,'SCAL14',scal14,0.4D0)
 131       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 132       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 133       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 134       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 135       call reada(weightcard,"D0CM",d0cm,3.78d0)
 136       call reada(weightcard,"AKCM",akcm,15.1d0)
 137       call reada(weightcard,"AKTH",akth,11.0d0)
 138       call reada(weightcard,"AKCT",akct,12.0d0)
 139       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 140       call reada(weightcard,"V2SS",v2ss,7.61d0)
 141       call reada(weightcard,"V3SS",v3ss,13.7d0)
 142       call reada(weightcard,"EBR",ebr,-5.50D0)
 143       if (index(weightcard,'SOFT').gt.0) ipot=6
 144 C 12/1/95 Added weight for the multi-body term WCORR
 145       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 146       do i=1,maxres-1
 147         do j=i+1,maxres
 148           dyn_ssbond_ij(i,j)=1.0d300
 149         enddo
 150       enddo
 151       call reada(weightcard,"HT",Ht,0.0D0)
 152       if (dyn_ss) then
 153         ss_depth=ebr/wsc-0.25*eps(1,1)
 154         Ht=Ht/wsc-0.25*eps(1,1)
 155         akcm=akcm*wstrain/wsc
 156         akth=akth*wstrain/wsc
 157         akct=akct*wstrain/wsc
 158         v1ss=v1ss*wstrain/wsc
 159         v2ss=v2ss*wstrain/wsc
 160         v3ss=v3ss*wstrain/wsc
 161       else
 162         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 163       endif
 164       write (iout,'(/a)') "Disulfide bridge parameters:"
 165       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 166         write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
 167       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 168       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 169       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 170      & ' v3ss:',v3ss
 171       write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
 172       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 173       weights(1)=wsc
 174       weights(2)=wscp
 175       weights(3)=welec
 176       weights(4)=wcorr
 177       weights(5)=wcorr5
 178       weights(6)=wcorr6
 179       weights(7)=wel_loc
 180       weights(8)=wturn3
 181       weights(9)=wturn4
 182       weights(10)=wturn6
 183       weights(11)=wang
 184       weights(12)=wscloc
 185       weights(13)=wtor
 186       weights(14)=wtor_d
 187       weights(15)=wstrain
 188       weights(16)=wvdwpp
 189       weights(17)=wbond
 190       weights(18)=scal14
 191       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 192      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 193      &  wturn4,wturn6,wsccor
 194    10 format (/'Energy-term weights (unscaled):'//
 195      & 'WSCC=   ',f10.6,' (SC-SC)'/
 196      & 'WSCP=   ',f10.6,' (SC-p)'/
 197      & 'WELEC=  ',f10.6,' (p-p electr)'/
 198      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 199      & 'WBOND=  ',f10.6,' (stretching)'/
 200      & 'WANG=   ',f10.6,' (bending)'/
 201      & 'WSCLOC= ',f10.6,' (SC local)'/
 202      & 'WTOR=   ',f10.6,' (torsional)'/
 203      & 'WTORD=  ',f10.6,' (double torsional)'/
 204      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 205      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 206      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 207      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 208      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 209      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 210      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 211      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 212      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 213       if (wcorr4.gt.0.0d0) then
 214         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 215      &   'between contact pairs of peptide groups'
 216         write (iout,'(2(a,f5.3/))')
 217      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 218      &  'Range of quenching the correlation terms:',2*delt_corr
 219       else if (wcorr.gt.0.0d0) then
 220         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 221      &   'between contact pairs of peptide groups'
 222       endif
 223       write (iout,'(a,f8.3)')
 224      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 225       write (iout,'(a,f8.3)')
 226      &  'General scaling factor of SC-p interactions:',scalscp
 227       r0_corr=cutoff_corr-delt_corr
 228       do i=1,20
 229         aad(i,1)=scalscp*aad(i,1)
 230         aad(i,2)=scalscp*aad(i,2)
 231         bad(i,1)=scalscp*bad(i,1)
 232         bad(i,2)=scalscp*bad(i,2)
 233       enddo
 234
 235       call flush(iout)
 236       print *,'indpdb=',indpdb,' pdbref=',pdbref
 237
 238 C Read sequence if not taken from the pdb file.
 239       if (iscode.gt.0) then
 240         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 241       else
 242         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 243       endif
 244 C Convert sequence to numeric code
 245       do i=1,nres
 246         itype(i)=rescode(i,sequence(i),iscode)
 247       enddo
 248       print *,nres
 249       print '(20i4)',(itype(i),i=1,nres)
 250
 251       do i=1,nres
 252 #ifdef PROCOR
 253         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 254 #else
 255         if (itype(i).eq.21) then
 256 #endif
 257           itel(i)=0
 258 #ifdef PROCOR
 259         else if (itype(i+1).ne.20) then
 260 #else
 261         else if (itype(i).ne.20) then
 262 #endif
 263           itel(i)=1
 264         else
 265           itel(i)=2
 266         endif
 267       enddo
 268       write (iout,*) "ITEL"
 269       do i=1,nres-1
 270         write (iout,*) i,itype(i),itel(i)
 271       enddo
 272
 273       print *,'Call Read_Bridge.'
 274       call read_bridge
 275       if (with_dihed_constr) then
 276
 277       read (inp,*) ndih_constr
 278       if (ndih_constr.gt.0) then
 279         read (inp,*) ftors
 280         write (iout,*) 'FTORS',ftors
 281         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 282         write (iout,*)
 283      &   'There are',ndih_constr,' constraints on phi angles.'
 284         do i=1,ndih_constr
 285           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 286         enddo
 287         do i=1,ndih_constr
 288           phi0(i)=deg2rad*phi0(i)
 289           drange(i)=deg2rad*drange(i)
 290         enddo
 291       endif
 292
 293       endif
 294
 295       nnt=1
 296       nct=nres
 297       print *,'NNT=',NNT,' NCT=',NCT
 298       if (itype(1).eq.21) nnt=2
 299       if (itype(nres).eq.21) nct=nct-1
 300       if (nstart.lt.nnt) nstart=nnt
 301       if (nend.gt.nct .or. nend.eq.0) nend=nct
 302       write (iout,*) "nstart",nstart," nend",nend
 303       nres0=nres
 304
 305       if (constr_homology.gt.0) then
 306         call read_constr_homology
 307       endif
 308
 309 c      if (pdbref) then
 310 c        read(inp,'(a)') pdbfile
 311 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 312 c        open(ipdbin,file=pdbfile,status='old',err=33)
 313 c        goto 34
 314 c  33    write (iout,'(a)') 'Error opening PDB file.'
 315 c        stop
 316 c  34    continue
 317 c        print *,'Begin reading pdb data'
 318 c        call readpdb
 319 c        print *,'Finished reading pdb data'
 320 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 321 c        do i=1,nres
 322 c          itype_pdb(i)=itype(i)
 323 c        enddo
 324 c        close (ipdbin)
 325 c        write (iout,'(a,i3)') 'nsup=',nsup
 326 c        nstart_seq=nnt
 327 c        if (nsup.le.(nct-nnt+1)) then
 328 c          do i=0,nct-nnt+1-nsup
 329 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 330 c              nstart_seq=nnt+i
 331 c              goto 111
 332 c            endif
 333 c          enddo
 334 c          write (iout,'(a)')
 335 c     &            'Error - sequences to be superposed do not match.'
 336 c          stop
 337 c        else
 338 c          do i=0,nsup-(nct-nnt+1)
 339 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 340 c     &      then
 341 c              nstart_sup=nstart_sup+i
 342 c              nsup=nct-nnt+1
 343 c              goto 111
 344 c            endif
 345 c          enddo
 346 c          write (iout,'(a)')
 347 c     &            'Error - sequences to be superposed do not match.'
 348 c        endif
 349 c  111   continue
 350 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 351 c     &                 ' nstart_seq=',nstart_seq
 352 c      endif
 353       call init_int_table
 354       call setup_var
 355       write (iout,*) "molread: REFSTR",refstr
 356       if (refstr) then
 357         if (.not.pdbref) then
 358           call read_angles(inp,*38)
 359           goto 39
 360    38     write (iout,'(a)') 'Error reading reference structure.'
 361 #ifdef MPL
 362           call mp_stopall(Error_Msg)
 363 #else
 364           stop 'Error reading reference structure'
 365 #endif
 366    39     call chainbuild
 367           nstart_sup=nnt
 368           nstart_seq=nnt
 369           nsup=nct-nnt+1
 370           do i=1,2*nres
 371             do j=1,3
 372               cref(j,i)=c(j,i)
 373             enddo
 374           enddo
 375         endif
 376         call contact(.true.,ncont_ref,icont_ref)
 377       endif
 378 c Read distance restraints
 379       if (constr_dist.gt.0) then
 380         call read_dist_constr
 381         call hpb_partition
 382       endif
 383       return
 384       end
 385 c-----------------------------------------------------------------------------
 386       logical function seq_comp(itypea,itypeb,length)
 387       implicit none
 388       integer length,itypea(length),itypeb(length)
 389       integer i
 390       do i=1,length
 391         if (itypea(i).ne.itypeb(i)) then
 392           seq_comp=.false.
 393           return
 394         endif
 395       enddo
 396       seq_comp=.true.
 397       return
 398       end
 399 c-----------------------------------------------------------------------------
 400       subroutine read_bridge
 401 C Read information about disulfide bridges.
 402       implicit none
 403       include 'DIMENSIONS'
 404       include 'COMMON.IOUNITS'
 405       include 'COMMON.GEO'
 406       include 'COMMON.VAR'
 407       include 'COMMON.INTERACT'
 408       include 'COMMON.LOCAL'
 409       include 'COMMON.NAMES'
 410       include 'COMMON.CHAIN'
 411       include 'COMMON.FFIELD'
 412       include 'COMMON.SBRIDGE'
 413       include 'COMMON.HEADER'
 414       include 'COMMON.CONTROL'
 415       include 'COMMON.TIME1'
 416 #ifdef MPL
 417       include 'COMMON.INFO'
 418 #endif
 419       integer i,j
 420 C Read bridging residues.
 421       read (inp,*) ns,(iss(i),i=1,ns)
 422       write(iout,*)'ns=',ns
 423 C Check whether the specified bridging residues are cystines.
 424       do i=1,ns
 425         if (itype(iss(i)).ne.1) then
 426           write (iout,'(2a,i3,a)')
 427      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 428      &   ' can form a disulfide bridge?!!!'
 429           write (*,'(2a,i3,a)')
 430      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 431      &   ' can form a disulfide bridge?!!!'
 432 #ifdef MPL
 433          call mp_stopall(error_msg)
 434 #else
 435          stop
 436 #endif
 437         endif
 438       enddo
 439 C Read preformed bridges.
 440       if (ns.gt.0) then
 441       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 442       if (nss.gt.0) then
 443         nhpb=nss
 444 C Check if the residues involved in bridges are in the specified list of
 445 C bridging residues.
 446         do i=1,nss
 447           do j=1,i-1
 448             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 449      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 450               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 451      &      ' contains residues present in other pairs.'
 452               write (*,'(a,i3,a)') 'Disulfide pair',i,
 453      &      ' contains residues present in other pairs.'
 454 #ifdef MPL
 455               call mp_stopall(error_msg)
 456 #else
 457               stop
 458 #endif
 459             endif
 460           enddo
 461           do j=1,ns
 462             if (ihpb(i).eq.iss(j)) goto 10
 463           enddo
 464           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 465    10     continue
 466           do j=1,ns
 467             if (jhpb(i).eq.iss(j)) goto 20
 468           enddo
 469           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 470    20     continue
 471 c          dhpb(i)=dbr
 472 c          forcon(i)=fbr
 473         enddo
 474         do i=1,nss
 475           ihpb(i)=ihpb(i)+nres
 476           jhpb(i)=jhpb(i)+nres
 477         enddo
 478       endif
 479       endif
 480       if (ns.gt.0.and.dyn_ss) then
 481           do i=nss+1,nhpb
 482             ihpb(i-nss)=ihpb(i)
 483             jhpb(i-nss)=jhpb(i)
 484             forcon(i-nss)=forcon(i)
 485             dhpb(i-nss)=dhpb(i)
 486           enddo
 487           nhpb=nhpb-nss
 488           nss=0
 489           call hpb_partition
 490           do i=1,ns
 491             dyn_ss_mask(iss(i))=.true.
 492 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 493           enddo
 494       endif
 495       print *, "Leaving brigde read"
 496       return
 497       end
 498 c----------------------------------------------------------------------------
 499       subroutine read_angles(kanal,*)
 500       implicit none
 501       include 'DIMENSIONS'
 502       include 'COMMON.GEO'
 503       include 'COMMON.VAR'
 504       include 'COMMON.CHAIN'
 505       include 'COMMON.IOUNITS'
 506       integer i,kanal
 507       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 508       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 509       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 510       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 511       do i=1,nres
 512         theta(i)=deg2rad*theta(i)
 513         phi(i)=deg2rad*phi(i)
 514         alph(i)=deg2rad*alph(i)
 515         omeg(i)=deg2rad*omeg(i)
 516       enddo
 517       return
 518    10 return1
 519       end
 520 c----------------------------------------------------------------------------
 521       subroutine reada(rekord,lancuch,wartosc,default)
 522       implicit none
 523       character*(*) rekord,lancuch
 524       double precision wartosc,default
 525       integer ilen,iread
 526       external ilen
 527       iread=index(rekord,lancuch)
 528       if (iread.eq.0) then
 529         wartosc=default
 530         return
 531       endif
 532       iread=iread+ilen(lancuch)+1
 533       read (rekord(iread:),*) wartosc
 534       return
 535       end
 536 c----------------------------------------------------------------------------
 537       subroutine multreada(rekord,lancuch,tablica,dim,default)
 538       implicit none
 539       integer dim,i
 540       double precision tablica(dim),default
 541       character*(*) rekord,lancuch
 542       integer ilen,iread
 543       external ilen
 544       do i=1,dim
 545         tablica(i)=default
 546       enddo
 547       iread=index(rekord,lancuch)
 548       if (iread.eq.0) return
 549       iread=iread+ilen(lancuch)+1
 550       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 551    10 return
 552       end
 553 c----------------------------------------------------------------------------
 554       subroutine readi(rekord,lancuch,wartosc,default)
 555       implicit none
 556       character*(*) rekord,lancuch
 557       integer wartosc,default
 558       integer ilen,iread
 559       external ilen
 560       iread=index(rekord,lancuch)
 561       if (iread.eq.0) then
 562         wartosc=default
 563         return
 564       endif
 565       iread=iread+ilen(lancuch)+1
 566       read (rekord(iread:),*) wartosc
 567       return
 568       end
 569 c----------------------------------------------------------------------------
 570       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 571       implicit none
 572       integer dim,i
 573       integer tablica(dim),default
 574       character*(*) rekord,lancuch
 575       character*80 aux
 576       integer ilen,iread
 577       external ilen
 578       do i=1,dim
 579         tablica(i)=default
 580       enddo
 581       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 582       if (iread.eq.0) return
 583       iread=iread+ilen(lancuch)+1
 584       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 585    10 return
 586       end
 587 c----------------------------------------------------------------------------
 588       subroutine card_concat(card)
 589       include 'DIMENSIONS'
 590       include 'COMMON.IOUNITS'
 591       character*(*) card
 592       character*80 karta,ucase
 593       external ilen
 594       read (inp,'(a)') karta
 595       karta=ucase(karta)
 596       card=' '
 597       do while (karta(80:80).eq.'&')
 598         card=card(:ilen(card)+1)//karta(:79)
 599         read (inp,'(a)') karta
 600         karta=ucase(karta)
 601       enddo
 602       card=card(:ilen(card)+1)//karta
 603       return
 604       end
 605 c----------------------------------------------------------------------------
 606       subroutine openunits
 607       implicit none
 608       include 'DIMENSIONS'
 609 #ifdef MPI
 610       include "mpif.h"
 611       character*3 liczba
 612       include "COMMON.MPI"
 613 #endif
 614       include 'COMMON.IOUNITS'
 615       include 'COMMON.CONTROL'
 616       integer lenpre,lenpot,ilen
 617       external ilen
 618       character*16 cformat,cprint
 619       character*16 ucase
 620       integer lenint,lenout
 621       call getenv('INPUT',prefix)
 622       call getenv('OUTPUT',prefout)
 623       call getenv('INTIN',prefintin)
 624       call getenv('COORD',cformat)
 625       call getenv('PRINTCOOR',cprint)
 626       call getenv('SCRATCHDIR',scratchdir)
 627       from_bx=.true.
 628       from_cx=.false.
 629       if (index(ucase(cformat),'CX').gt.0) then
 630         from_cx=.true.
 631         from_bx=.false.
 632       endif
 633       from_cart=.true.
 634       lenpre=ilen(prefix)
 635       lenout=ilen(prefout)
 636       lenint=ilen(prefintin)
 637 C Get the names and open the input files
 638       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 639 #ifdef MPI
 640       write (liczba,'(bz,i3.3)') me
 641       outname=prefout(:lenout)//'_clust.out_'//liczba
 642 #else
 643       outname=prefout(:lenout)//'_clust.out'
 644 #endif
 645       if (from_bx) then
 646         intinname=prefintin(:lenint)//'.bx'
 647       else if (from_cx) then
 648         intinname=prefintin(:lenint)//'.cx'
 649       else
 650         intinname=prefintin(:lenint)//'.int'
 651       endif
 652       rmsname=prefintin(:lenint)//'.rms'
 653       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 654      &       status='unknown')
 655       open (jrms,file=rmsname,status='unknown')
 656       open(iout,file=outname,status='unknown')
 657 C Get parameter filenames and open the parameter files.
 658       call getenv('BONDPAR',bondname)
 659       open (ibond,file=bondname,status='old')
 660       call getenv('THETPAR',thetname)
 661       open (ithep,file=thetname,status='old')
 662       call getenv('ROTPAR',rotname)
 663       open (irotam,file=rotname,status='old')
 664       call getenv('TORPAR',torname)
 665       open (itorp,file=torname,status='old')
 666       call getenv('TORDPAR',tordname)
 667       open (itordp,file=tordname,status='old')
 668       call getenv('FOURIER',fouriername)
 669       open (ifourier,file=fouriername,status='old')
 670       call getenv('ELEPAR',elename)
 671       open (ielep,file=elename,status='old')
 672       call getenv('SIDEPAR',sidename)
 673       open (isidep,file=sidename,status='old')
 674       call getenv('SIDEP',sidepname)
 675       open (isidep1,file=sidepname,status="old")
 676       call getenv('SCCORPAR',sccorname)
 677       open (isccor,file=sccorname,status="old")
 678 #ifndef OLDSCP
 679 C
 680 C 8/9/01 In the newest version SCp interaction constants are read from a file
 681 C Use -DOLDSCP to use hard-coded constants instead.
 682 C
 683       call getenv('SCPPAR',scpname)
 684       open (iscpp,file=scpname,status='old')
 685 #endif
 686       return
 687       end
 688 c-------------------------------------------------------------------------------
 689       subroutine read_dist_constr
 690       implicit real*8 (a-h,o-z)
 691       include 'DIMENSIONS'
 692       include 'COMMON.CONTROL'
 693       include 'COMMON.CHAIN'
 694       include 'COMMON.IOUNITS'
 695       include 'COMMON.SBRIDGE'
 696       integer ifrag_(2,100),ipair_(2,100)
 697       double precision wfrag_(100),wpair_(100)
 698       character*500 controlcard
 699 c      write (iout,*) "Calling read_dist_constr"
 700 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 701       call card_concat(controlcard)
 702 c      call flush(iout)
 703 c      write (iout,'(a)') controlcard
 704       call readi(controlcard,"NFRAG",nfrag_,0)
 705       call readi(controlcard,"NPAIR",npair_,0)
 706       call readi(controlcard,"NDIST",ndist_,0)
 707       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 708       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 709       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 710       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 711       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 712       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 713       write (iout,*) "IFRAG"
 714       do i=1,nfrag_
 715         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 716       enddo
 717       write (iout,*) "IPAIR"
 718       do i=1,npair_
 719         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 720       enddo
 721       call flush(iout)
 722       if (.not.refstr .and. nfrag_.gt.0) then
 723         write (iout,*)
 724      &  "ERROR: no reference structure to compute distance restraints"
 725         write (iout,*)
 726      &  "Restraints must be specified explicitly (NDIST=number)"
 727         stop
 728       endif
 729       if (nfrag_.lt.2 .and. npair_.gt.0) then
 730         write (iout,*) "ERROR: Less than 2 fragments specified",
 731      &   " but distance restraints between pairs requested"
 732         stop
 733       endif
 734       call flush(iout)
 735       do i=1,nfrag_
 736         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 737         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 738      &    ifrag_(2,i)=nstart_sup+nsup-1
 739 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 740         call flush(iout)
 741         if (wfrag_(i).gt.0.0d0) then
 742         do j=ifrag_(1,i),ifrag_(2,i)-1
 743           do k=j+1,ifrag_(2,i)
 744             write (iout,*) "j",j," k",k
 745             ddjk=dist(j,k)
 746             if (constr_dist.eq.1) then
 747             nhpb=nhpb+1
 748             ihpb(nhpb)=j
 749             jhpb(nhpb)=k
 750               dhpb(nhpb)=ddjk
 751             forcon(nhpb)=wfrag_(i)
 752             else if (constr_dist.eq.2) then
 753               if (ddjk.le.dist_cut) then
 754                 nhpb=nhpb+1
 755                 ihpb(nhpb)=j
 756                 jhpb(nhpb)=k
 757                 dhpb(nhpb)=ddjk
 758                 forcon(nhpb)=wfrag_(i)
 759               endif
 760             else
 761               nhpb=nhpb+1
 762               ihpb(nhpb)=j
 763               jhpb(nhpb)=k
 764               dhpb(nhpb)=ddjk
 765               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 766             endif
 767             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 768      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 769           enddo
 770         enddo
 771         endif
 772       enddo
 773       do i=1,npair_
 774         if (wpair_(i).gt.0.0d0) then
 775         ii = ipair_(1,i)
 776         jj = ipair_(2,i)
 777         if (ii.gt.jj) then
 778           itemp=ii
 779           ii=jj
 780           jj=itemp
 781         endif
 782         do j=ifrag_(1,ii),ifrag_(2,ii)
 783           do k=ifrag_(1,jj),ifrag_(2,jj)
 784             nhpb=nhpb+1
 785             ihpb(nhpb)=j
 786             jhpb(nhpb)=k
 787             forcon(nhpb)=wpair_(i)
 788             dhpb(nhpb)=dist(j,k)
 789             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 790      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 791           enddo
 792         enddo
 793         endif
 794       enddo
 795       do i=1,ndist_
 796         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 797      &     ibecarb(i),forcon(nhpb+1)
 798         if (forcon(nhpb+1).gt.0.0d0) then
 799           nhpb=nhpb+1
 800           if (ibecarb(i).gt.0) then
 801             ihpb(i)=ihpb(i)+nres
 802             jhpb(i)=jhpb(i)+nres
 803           endif
 804           if (dhpb(nhpb).eq.0.0d0)
 805      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 806         endif
 807       enddo
 808       do i=1,nhpb
 809           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 810      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 811       enddo
 812       call flush(iout)
 813       return
 814       end
 815
 816 c====-------------------------------------------------------------------
 817       subroutine read_constr_homology
 818
 819       include 'DIMENSIONS'
 820 #ifdef MPI
 821       include 'mpif.h'
 822 #endif
 823       include 'COMMON.CONTROL'
 824       include 'COMMON.CHAIN'
 825       include 'COMMON.IOUNITS'
 826       include 'COMMON.INTERACT'
 827       include 'COMMON.GEO'
 828       double precision odl_temp,sigma_odl_temp
 829       common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 830      &    sigma_odl_temp(maxres,maxres,max_template)
 831       character*2 kic2
 832       character*24 model_ki_dist, model_ki_angle
 833       character*500 controlcard
 834       integer ki, i, j, k, l
 835       logical lprn /.true./
 836
 837       call card_concat(controlcard)
 838       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 839       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 840       write (iout,*) "wga_dist",waga_dist," waga_angle",waga_angle
 841
 842       lim_odl=0
 843       lim_dih=0
 844       do i=1,nres
 845         do j=i+2,nres
 846           do ki=1,constr_homology
 847             sigma_odl_temp(i,j,ki)=0.0d0
 848             odl_temp(i,j,ki)=0.0d0
 849           enddo
 850         enddo
 851       enddo
 852       do i=1,nres-3
 853         do ki=1,constr_homology
 854           dih(ki,i)=0.0d0
 855           sigma_dih(ki,i)=0.0d0
 856         enddo
 857       enddo
 858       do ki=1,constr_homology
 859           write(kic2,'(i2)') ki
 860           if (ki.le.9) kic2="0"//kic2(2:2)
 861
 862           model_ki_dist="model"//kic2//".dist"
 863           model_ki_angle="model"//kic2//".angle"
 864         open (ientin,file=model_ki_dist,status='old')
 865         do irec=1,maxdim !petla do czytania wiezow na odleglosc
 866           read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
 867      &       sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
 868           odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
 869           sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
 870      &     sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
 871         enddo
 872  1401 continue
 873         close (ientin)
 874         open (ientin,file=model_ki_angle,status='old')
 875         do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
 876           read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
 877      &      sigma_dih(ki,i+nnt-1)
 878           write (iout,*) i, j, k,l,dih(ki,i+nnt-1),
 879      &      sigma_dih(ki,i+nnt-1)
 880           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
 881           sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
 882         enddo
 883  1402 continue
 884         close (ientin)
 885       enddo
 886       ii=0
 887       write (iout,*) "nnt",nnt," nct",nct
 888       do i=nnt,nct-2
 889         do j=i+2,nct
 890           ki=1
 891 c          write (iout,*) "i",i," j",j," constr_homology",constr_homology
 892           do while (ki.le.constr_homology .and.
 893      &        sigma_odl_temp(i,j,ki).le.0.0d0)
 894 c            write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
 895             ki=ki+1
 896           enddo
 897 c          write (iout,*) "ki",ki
 898           if (ki.gt.constr_homology) cycle
 899           ii=ii+1
 900           ires_homo(ii)=i
 901           jres_homo(ii)=j
 902           do ki=1,constr_homology
 903             odl(ki,ii)=odl_temp(i,j,ki)
 904             sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
 905           enddo
 906         enddo
 907       enddo
 908       lim_odl=ii
 909       if (constr_homology.gt.0) call homology_partition
 910 c Print restraints
 911       if (.not.lprn) return
 912       write (iout,*) "Distance restraints from templates"
 913       do ii=1,lim_odl
 914         write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
 915      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
 916       enddo
 917       write (iout,*) "Dihedral angle restraints from templates"
 918       do i=nnt,lim_dih
 919         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
 920      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 921       enddo
 922 c      write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
 923 c      write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
 924
 925       return
 926       end
 927
 928