source/cluster/wham/src-restraint-MP/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       if (constr_homology.gt.0) then
 296         call read_constr_homology
 297       endif
 298
 299       nnt=1
 300       nct=nres
 301       print *,'NNT=',NNT,' NCT=',NCT
 302       if (itype(1).eq.21) nnt=2
 303       if (itype(nres).eq.21) nct=nct-1
 304       if (nstart.lt.nnt) nstart=nnt
 305       if (nend.gt.nct .or. nend.eq.0) nend=nct
 306       write (iout,*) "nstart",nstart," nend",nend
 307       nres0=nres
 308 c      if (pdbref) then
 309 c        read(inp,'(a)') pdbfile
 310 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 311 c        open(ipdbin,file=pdbfile,status='old',err=33)
 312 c        goto 34
 313 c  33    write (iout,'(a)') 'Error opening PDB file.'
 314 c        stop
 315 c  34    continue
 316 c        print *,'Begin reading pdb data'
 317 c        call readpdb
 318 c        print *,'Finished reading pdb data'
 319 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 320 c        do i=1,nres
 321 c          itype_pdb(i)=itype(i)
 322 c        enddo
 323 c        close (ipdbin)
 324 c        write (iout,'(a,i3)') 'nsup=',nsup
 325 c        nstart_seq=nnt
 326 c        if (nsup.le.(nct-nnt+1)) then
 327 c          do i=0,nct-nnt+1-nsup
 328 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 329 c              nstart_seq=nnt+i
 330 c              goto 111
 331 c            endif
 332 c          enddo
 333 c          write (iout,'(a)')
 334 c     &            'Error - sequences to be superposed do not match.'
 335 c          stop
 336 c        else
 337 c          do i=0,nsup-(nct-nnt+1)
 338 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 339 c     &      then
 340 c              nstart_sup=nstart_sup+i
 341 c              nsup=nct-nnt+1
 342 c              goto 111
 343 c            endif
 344 c          enddo
 345 c          write (iout,'(a)')
 346 c     &            'Error - sequences to be superposed do not match.'
 347 c        endif
 348 c  111   continue
 349 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 350 c     &                 ' nstart_seq=',nstart_seq
 351 c      endif
 352       call init_int_table
 353       call setup_var
 354       write (iout,*) "molread: REFSTR",refstr
 355       if (refstr) then
 356         if (.not.pdbref) then
 357           call read_angles(inp,*38)
 358           goto 39
 359    38     write (iout,'(a)') 'Error reading reference structure.'
 360 #ifdef MPL
 361           call mp_stopall(Error_Msg)
 362 #else
 363           stop 'Error reading reference structure'
 364 #endif
 365    39     call chainbuild
 366           nstart_sup=nnt
 367           nstart_seq=nnt
 368           nsup=nct-nnt+1
 369           do i=1,2*nres
 370             do j=1,3
 371               cref(j,i)=c(j,i)
 372             enddo
 373           enddo
 374         endif
 375         call contact(.true.,ncont_ref,icont_ref)
 376       endif
 377 c Read distance restraints
 378       if (constr_dist.gt.0) then
 379         call read_dist_constr
 380         call hpb_partition
 381       endif
 382       return
 383       end
 384 c-----------------------------------------------------------------------------
 385       logical function seq_comp(itypea,itypeb,length)
 386       implicit none
 387       integer length,itypea(length),itypeb(length)
 388       integer i
 389       do i=1,length
 390         if (itypea(i).ne.itypeb(i)) then
 391           seq_comp=.false.
 392           return
 393         endif
 394       enddo
 395       seq_comp=.true.
 396       return
 397       end
 398 c-----------------------------------------------------------------------------
 399       subroutine read_bridge
 400 C Read information about disulfide bridges.
 401       implicit none
 402       include 'DIMENSIONS'
 403       include 'COMMON.IOUNITS'
 404       include 'COMMON.GEO'
 405       include 'COMMON.VAR'
 406       include 'COMMON.INTERACT'
 407       include 'COMMON.LOCAL'
 408       include 'COMMON.NAMES'
 409       include 'COMMON.CHAIN'
 410       include 'COMMON.FFIELD'
 411       include 'COMMON.SBRIDGE'
 412       include 'COMMON.HEADER'
 413       include 'COMMON.CONTROL'
 414       include 'COMMON.TIME1'
 415 #ifdef MPL
 416       include 'COMMON.INFO'
 417 #endif
 418       integer i,j
 419 C Read bridging residues.
 420       read (inp,*) ns,(iss(i),i=1,ns)
 421       print *,'ns=',ns
 422 C Check whether the specified bridging residues are cystines.
 423       do i=1,ns
 424         if (itype(iss(i)).ne.1) then
 425           write (iout,'(2a,i3,a)')
 426      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 427      &   ' can form a disulfide bridge?!!!'
 428           write (*,'(2a,i3,a)')
 429      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 430      &   ' can form a disulfide bridge?!!!'
 431 #ifdef MPL
 432          call mp_stopall(error_msg)
 433 #else
 434          stop
 435 #endif
 436         endif
 437       enddo
 438 C Read preformed bridges.
 439       if (ns.gt.0) then
 440       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 441       if (nss.gt.0) then
 442         nhpb=nss
 443 C Check if the residues involved in bridges are in the specified list of
 444 C bridging residues.
 445         do i=1,nss
 446           do j=1,i-1
 447             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 448      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 449               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 450      &      ' contains residues present in other pairs.'
 451               write (*,'(a,i3,a)') 'Disulfide pair',i,
 452      &      ' contains residues present in other pairs.'
 453 #ifdef MPL
 454               call mp_stopall(error_msg)
 455 #else
 456               stop
 457 #endif
 458             endif
 459           enddo
 460           do j=1,ns
 461             if (ihpb(i).eq.iss(j)) goto 10
 462           enddo
 463           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 464    10     continue
 465           do j=1,ns
 466             if (jhpb(i).eq.iss(j)) goto 20
 467           enddo
 468           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 469    20     continue
 470 c          dhpb(i)=dbr
 471 c          forcon(i)=fbr
 472         enddo
 473         do i=1,nss
 474           ihpb(i)=ihpb(i)+nres
 475           jhpb(i)=jhpb(i)+nres
 476         enddo
 477       endif
 478       endif
 479       if (ns.gt.0.and.dyn_ss) then
 480           do i=nss+1,nhpb
 481             ihpb(i-nss)=ihpb(i)
 482             jhpb(i-nss)=jhpb(i)
 483             forcon(i-nss)=forcon(i)
 484             dhpb(i-nss)=dhpb(i)
 485           enddo
 486           nhpb=nhpb-nss
 487           nss=0
 488           call hpb_partition
 489           do i=1,ns
 490             dyn_ss_mask(iss(i))=.true.
 491 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 492           enddo
 493       endif
 494       print *, "Leaving brigde read"
 495       return
 496       end
 497 c----------------------------------------------------------------------------
 498       subroutine read_angles(kanal,*)
 499       implicit none
 500       include 'DIMENSIONS'
 501       include 'COMMON.GEO'
 502       include 'COMMON.VAR'
 503       include 'COMMON.CHAIN'
 504       include 'COMMON.IOUNITS'
 505       integer i,kanal
 506       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 507       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 508       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 509       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 510       do i=1,nres
 511         theta(i)=deg2rad*theta(i)
 512         phi(i)=deg2rad*phi(i)
 513         alph(i)=deg2rad*alph(i)
 514         omeg(i)=deg2rad*omeg(i)
 515       enddo
 516       return
 517    10 return1
 518       end
 519 c----------------------------------------------------------------------------
 520       subroutine reada(rekord,lancuch,wartosc,default)
 521       implicit none
 522       character*(*) rekord,lancuch
 523       double precision wartosc,default
 524       integer ilen,iread
 525       external ilen
 526       iread=index(rekord,lancuch)
 527       if (iread.eq.0) then
 528         wartosc=default
 529         return
 530       endif
 531       iread=iread+ilen(lancuch)+1
 532       read (rekord(iread:),*) wartosc
 533       return
 534       end
 535 c----------------------------------------------------------------------------
 536       subroutine multreada(rekord,lancuch,tablica,dim,default)
 537       implicit none
 538       integer dim,i
 539       double precision tablica(dim),default
 540       character*(*) rekord,lancuch
 541       integer ilen,iread
 542       external ilen
 543       do i=1,dim
 544         tablica(i)=default
 545       enddo
 546       iread=index(rekord,lancuch)
 547       if (iread.eq.0) return
 548       iread=iread+ilen(lancuch)+1
 549       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 550    10 return
 551       end
 552 c----------------------------------------------------------------------------
 553       subroutine readi(rekord,lancuch,wartosc,default)
 554       implicit none
 555       character*(*) rekord,lancuch
 556       integer wartosc,default
 557       integer ilen,iread
 558       external ilen
 559       iread=index(rekord,lancuch)
 560       if (iread.eq.0) then
 561         wartosc=default
 562         return
 563       endif
 564       iread=iread+ilen(lancuch)+1
 565       read (rekord(iread:),*) wartosc
 566       return
 567       end
 568 c----------------------------------------------------------------------------
 569       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 570       implicit none
 571       integer dim,i
 572       integer tablica(dim),default
 573       character*(*) rekord,lancuch
 574       character*80 aux
 575       integer ilen,iread
 576       external ilen
 577       do i=1,dim
 578         tablica(i)=default
 579       enddo
 580       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 581       if (iread.eq.0) return
 582       iread=iread+ilen(lancuch)+1
 583       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 584    10 return
 585       end
 586 c----------------------------------------------------------------------------
 587       subroutine card_concat(card)
 588       include 'DIMENSIONS'
 589       include 'COMMON.IOUNITS'
 590       character*(*) card
 591       character*80 karta,ucase
 592       external ilen
 593       read (inp,'(a)') karta
 594       karta=ucase(karta)
 595       card=' '
 596       do while (karta(80:80).eq.'&')
 597         card=card(:ilen(card)+1)//karta(:79)
 598         read (inp,'(a)') karta
 599         karta=ucase(karta)
 600       enddo
 601       card=card(:ilen(card)+1)//karta
 602       return
 603       end
 604 c----------------------------------------------------------------------------
 605       subroutine openunits
 606       implicit none
 607       include 'DIMENSIONS'
 608 #ifdef MPI
 609       include "mpif.h"
 610       character*3 liczba
 611       include "COMMON.MPI"
 612 #endif
 613       include 'COMMON.IOUNITS'
 614       include 'COMMON.CONTROL'
 615       integer lenpre,lenpot,ilen
 616       external ilen
 617       character*16 cformat,cprint
 618       character*16 ucase
 619       integer lenint,lenout
 620       call getenv('INPUT',prefix)
 621       call getenv('OUTPUT',prefout)
 622       call getenv('INTIN',prefintin)
 623       call getenv('COORD',cformat)
 624       call getenv('PRINTCOOR',cprint)
 625       call getenv('SCRATCHDIR',scratchdir)
 626       from_bx=.true.
 627       from_cx=.false.
 628       if (index(ucase(cformat),'CX').gt.0) then
 629         from_cx=.true.
 630         from_bx=.false.
 631       endif
 632       from_cart=.true.
 633       lenpre=ilen(prefix)
 634       lenout=ilen(prefout)
 635       lenint=ilen(prefintin)
 636 C Get the names and open the input files
 637       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 638 #ifdef MPI
 639       write (liczba,'(bz,i3.3)') me
 640       outname=prefout(:lenout)//'_clust.out_'//liczba
 641 #else
 642       outname=prefout(:lenout)//'_clust.out'
 643 #endif
 644       if (from_bx) then
 645         intinname=prefintin(:lenint)//'.bx'
 646       else if (from_cx) then
 647         intinname=prefintin(:lenint)//'.cx'
 648       else
 649         intinname=prefintin(:lenint)//'.int'
 650       endif
 651       rmsname=prefintin(:lenint)//'.rms'
 652       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 653      &       status='unknown')
 654       open (jrms,file=rmsname,status='unknown')
 655       open(iout,file=outname,status='unknown')
 656 C Get parameter filenames and open the parameter files.
 657       call getenv('BONDPAR',bondname)
 658       open (ibond,file=bondname,status='old')
 659       call getenv('THETPAR',thetname)
 660       open (ithep,file=thetname,status='old')
 661       call getenv('ROTPAR',rotname)
 662       open (irotam,file=rotname,status='old')
 663       call getenv('TORPAR',torname)
 664       open (itorp,file=torname,status='old')
 665       call getenv('TORDPAR',tordname)
 666       open (itordp,file=tordname,status='old')
 667       call getenv('FOURIER',fouriername)
 668       open (ifourier,file=fouriername,status='old')
 669       call getenv('ELEPAR',elename)
 670       open (ielep,file=elename,status='old')
 671       call getenv('SIDEPAR',sidename)
 672       open (isidep,file=sidename,status='old')
 673       call getenv('SIDEP',sidepname)
 674       open (isidep1,file=sidepname,status="old")
 675       call getenv('SCCORPAR',sccorname)
 676       open (isccor,file=sccorname,status="old")
 677 #ifndef OLDSCP
 678 C
 679 C 8/9/01 In the newest version SCp interaction constants are read from a file
 680 C Use -DOLDSCP to use hard-coded constants instead.
 681 C
 682       call getenv('SCPPAR',scpname)
 683       open (iscpp,file=scpname,status='old')
 684 #endif
 685       return
 686       end
 687 c-------------------------------------------------------------------------------
 688       subroutine read_dist_constr
 689       implicit real*8 (a-h,o-z)
 690       include 'DIMENSIONS'
 691       include 'COMMON.CONTROL'
 692       include 'COMMON.CHAIN'
 693       include 'COMMON.IOUNITS'
 694       include 'COMMON.SBRIDGE'
 695       integer ifrag_(2,100),ipair_(2,100)
 696       double precision wfrag_(100),wpair_(100)
 697       character*500 controlcard
 698 c      write (iout,*) "Calling read_dist_constr"
 699 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 700       call card_concat(controlcard)
 701 c      call flush(iout)
 702 c      write (iout,'(a)') controlcard
 703       call readi(controlcard,"NFRAG",nfrag_,0)
 704       call readi(controlcard,"NPAIR",npair_,0)
 705       call readi(controlcard,"NDIST",ndist_,0)
 706       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 707       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 708       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 709       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 710       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 711       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 712       write (iout,*) "IFRAG"
 713       do i=1,nfrag_
 714         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 715       enddo
 716       write (iout,*) "IPAIR"
 717       do i=1,npair_
 718         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 719       enddo
 720       call flush(iout)
 721       if (.not.refstr .and. nfrag_.gt.0) then
 722         write (iout,*)
 723      &  "ERROR: no reference structure to compute distance restraints"
 724         write (iout,*)
 725      &  "Restraints must be specified explicitly (NDIST=number)"
 726         stop
 727       endif
 728       if (nfrag_.lt.2 .and. npair_.gt.0) then
 729         write (iout,*) "ERROR: Less than 2 fragments specified",
 730      &   " but distance restraints between pairs requested"
 731         stop
 732       endif
 733       call flush(iout)
 734       do i=1,nfrag_
 735         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 736         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 737      &    ifrag_(2,i)=nstart_sup+nsup-1
 738 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 739         call flush(iout)
 740         if (wfrag_(i).gt.0.0d0) then
 741         do j=ifrag_(1,i),ifrag_(2,i)-1
 742           do k=j+1,ifrag_(2,i)
 743             write (iout,*) "j",j," k",k
 744             ddjk=dist(j,k)
 745             if (constr_dist.eq.1) then
 746             nhpb=nhpb+1
 747             ihpb(nhpb)=j
 748             jhpb(nhpb)=k
 749               dhpb(nhpb)=ddjk
 750             forcon(nhpb)=wfrag_(i)
 751             else if (constr_dist.eq.2) then
 752               if (ddjk.le.dist_cut) then
 753                 nhpb=nhpb+1
 754                 ihpb(nhpb)=j
 755                 jhpb(nhpb)=k
 756                 dhpb(nhpb)=ddjk
 757                 forcon(nhpb)=wfrag_(i)
 758               endif
 759             else
 760               nhpb=nhpb+1
 761               ihpb(nhpb)=j
 762               jhpb(nhpb)=k
 763               dhpb(nhpb)=ddjk
 764               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 765             endif
 766             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 767      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 768           enddo
 769         enddo
 770         endif
 771       enddo
 772       do i=1,npair_
 773         if (wpair_(i).gt.0.0d0) then
 774         ii = ipair_(1,i)
 775         jj = ipair_(2,i)
 776         if (ii.gt.jj) then
 777           itemp=ii
 778           ii=jj
 779           jj=itemp
 780         endif
 781         do j=ifrag_(1,ii),ifrag_(2,ii)
 782           do k=ifrag_(1,jj),ifrag_(2,jj)
 783             nhpb=nhpb+1
 784             ihpb(nhpb)=j
 785             jhpb(nhpb)=k
 786             forcon(nhpb)=wpair_(i)
 787             dhpb(nhpb)=dist(j,k)
 788             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 789      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 790           enddo
 791         enddo
 792         endif
 793       enddo
 794       do i=1,ndist_
 795         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 796      &     ibecarb(i),forcon(nhpb+1)
 797         if (forcon(nhpb+1).gt.0.0d0) then
 798           nhpb=nhpb+1
 799           if (ibecarb(i).gt.0) then
 800             ihpb(i)=ihpb(i)+nres
 801             jhpb(i)=jhpb(i)+nres
 802           endif
 803           if (dhpb(nhpb).eq.0.0d0)
 804      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 805         endif
 806       enddo
 807       do i=1,nhpb
 808           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 809      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 810       enddo
 811       call flush(iout)
 812       return
 813       end
 814
 815 c====-------------------------------------------------------------------
 816
 817       subroutine read_constr_homology
 818
 819       include 'DIMENSIONS'
 820 #ifdef MPI
 821       include 'mpif.h'
 822 #endif
 823 c      include 'COMMON.SETUP'
 824       include 'COMMON.CONTROL'
 825       include 'COMMON.CHAIN'
 826       include 'COMMON.IOUNITS'
 827 c      include 'COMMON.MD'
 828
 829       character*2 kic2
 830       character*24 model_ki_dist, model_ki_angle
 831       character*500 controlcard
 832       integer ki, i, j, k, l
 833
 834 c      write(iout,*)"TEST1", waga_dist, waga_angle
 835 c      call card_concat(controlcard)
 836       call card_concat(controlcard)
 837 c      write(iout,*)"TEST1a", waga_dist, waga_angle
 838       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0)
 839 c      write(iout,*)"TEST1b", waga_dist, waga_angle
 840       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0)
 841
 842       write(iout,*)"TEST2", waga_dist, waga_angle
 843
 844       write(iout,*) "TEST_constr_homology_czytanie", constr_homology
 845       do ki=1,constr_homology
 846         if (constr_homology.ge.1) then
 847
 848           write(kic2,'(i2)') ki
 849 c          write(iout,*) "TEST KICA, HOMOL", kic2
 850           if (ki.le.9) kic2="0"//kic2(2:2)
 851 c          write(iout,*) "TEST KICA2, HOMOL", kic2
 852           model_ki_dist="model"//kic2//".dist"
 853           model_ki_angle="model"//kic2//".angle"
 854 c      write(iout,*) model_ki_dist, model_ki_angle
 855           open (1400+ki,file=model_ki_dist,status='old')
 856           open (1401+ki,file=model_ki_angle,status='old')
 857
 858         do irec=1,99999 !petla do czytania wiezow na odleglosc
 859           read (1400+ki,*,end=1401) i, j, odl(i,j,ki),sigma_odl(i,j,ki)
 860           lim_odl=i
 861         enddo
 862  1401 continue
 863 c      write(iout,*) "TEST_CZYTANIA__WIEZOW_ODL", i, j, odl(i,j,ki)
 864         do irec=1,99999 !petla do czytania wiezow na katach torsyjnych
 865           read (1401+ki,*,end=1402) i, j, k,l,dih(i,ki),sigma_dih(i,ki)
 866           lim_dih=i
 867 c          dih(i,ki)=dih(i,ki)
 868         enddo
 869  1402 continue
 870         endif
 871       enddo
 872
 873 c      write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
 874 c      write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
 875
 876
 877       return
 878       end
 879