source/analysis/src-M-prop/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       include 'COMMON.INTERACT'
  18       include "COMMON.SPLITELE"
  19       character*320 controlcard,ucase
  20 #ifdef MPL
  21       include 'COMMON.INFO'
  22 #endif
  23       integer i
  24
  25       read (INP,'(a80)') titel
  26       call card_concat(controlcard)
  27
  28       call readi(controlcard,'NRES',nres,0)
  29       call readi(controlcard,'RESCALE',rescale_mode,2)
  30       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
  31       write (iout,*) "DISTCHAINMAX",distchainmax
  32 C Reading the dimensions of box in x,y,z coordinates
  33       call reada(controlcard,'BOXX',boxxsize,100.0d0)
  34       call reada(controlcard,'BOXY',boxysize,100.0d0)
  35       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
  36 c Cutoff range for interactions
  37       call reada(controlcard,"R_CUT",r_cut,15.0d0)
  38       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
  39       refstr=(index(controlcard,'REFSTR').gt.0)
  40       write (iout,*) "REFSTR",refstr
  41       pdbref=(index(controlcard,'PDBREF').gt.0)
  42       iscode=index(controlcard,'ONE_LETTER')
  43       call readi(controlcard,'SYM',symetr,1)
  44       write (iout,*) 'sym', symetr
  45       call readi(controlcard,'NSTART',nstart,0)
  46       call readi(controlcard,'NEND',nend,0)
  47       caonly=(index(controlcard,'CA_ONLY').gt.0)
  48       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
  49       lside = index(controlcard,"SIDE").gt.0
  50       return
  51       end
  52 c--------------------------------------------------------------------------
  53       subroutine molread
  54 C
  55 C Read molecular data.
  56 C
  57       implicit none
  58       include 'DIMENSIONS'
  59       include 'COMMON.IOUNITS'
  60       include 'COMMON.GEO'
  61       include 'COMMON.VAR'
  62       include 'COMMON.INTERACT'
  63       include 'COMMON.LOCAL'
  64       include 'COMMON.NAMES'
  65       include 'COMMON.CHAIN'
  66       include 'COMMON.FFIELD'
  67       include 'COMMON.SBRIDGE'
  68       include 'COMMON.HEADER'
  69       include 'COMMON.CONTROL'
  70       include 'COMMON.CONTACTS'
  71       include 'COMMON.TIME1'
  72 #ifdef MPL
  73       include 'COMMON.INFO'
  74 #endif
  75       character*4 sequence(maxres)
  76       character*800 weightcard
  77       integer rescode
  78       double precision x(maxvar)
  79       integer itype_pdb(maxres)
  80       logical seq_comp
  81       integer i,j,kkk
  82 C
  83 C Body
  84 C
  85 C Read weights of the subsequent energy terms.
  86       call card_concat(weightcard)
  87       call reada(weightcard,'WSC',wsc,1.0d0)
  88       call reada(weightcard,'WLONG',wsc,wsc)
  89       call reada(weightcard,'WSCP',wscp,1.0d0)
  90       call reada(weightcard,'WELEC',welec,1.0D0)
  91       call reada(weightcard,'WVDWPP',wvdwpp,welec)
  92       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
  93       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
  94       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
  95       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
  96       call reada(weightcard,'WTURN3',wturn3,1.0D0)
  97       call reada(weightcard,'WTURN4',wturn4,1.0D0)
  98       call reada(weightcard,'WTURN6',wturn6,1.0D0)
  99       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 100       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 101       call reada(weightcard,'WBOND',wbond,1.0D0)
 102       call reada(weightcard,'WTOR',wtor,1.0D0)
 103       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 104       call reada(weightcard,'WANG',wang,1.0D0)
 105       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 106       call reada(weightcard,'SCAL14',scal14,0.4D0)
 107       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 108       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 109       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 110       if (index(weightcard,'SOFT').gt.0) ipot=6
 111 C 12/1/95 Added weight for the multi-body term WCORR
 112       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 113       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 114       weights(1)=wsc
 115       weights(2)=wscp
 116       weights(3)=welec
 117       weights(4)=wcorr
 118       weights(5)=wcorr5
 119       weights(6)=wcorr6
 120       weights(7)=wel_loc
 121       weights(8)=wturn3
 122       weights(9)=wturn4
 123       weights(10)=wturn6
 124       weights(11)=wang
 125       weights(12)=wscloc
 126       weights(13)=wtor
 127       weights(14)=wtor_d
 128       weights(15)=wstrain
 129       weights(16)=wvdwpp
 130       weights(17)=wbond
 131       weights(18)=scal14
 132       weights(19)=wsccor
 133       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 134      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 135      &  wturn4,wturn6,wsccor
 136    10 format (/'Energy-term weights (unscaled):'//
 137      & 'WSCC=   ',f10.6,' (SC-SC)'/
 138      & 'WSCP=   ',f10.6,' (SC-p)'/
 139      & 'WELEC=  ',f10.6,' (p-p electr)'/
 140      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 141      & 'WBOND=  ',f10.6,' (stretching)'/
 142      & 'WANG=   ',f10.6,' (bending)'/
 143      & 'WSCLOC= ',f10.6,' (SC local)'/
 144      & 'WTOR=   ',f10.6,' (torsional)'/
 145      & 'WTORD=  ',f10.6,' (double torsional)'/
 146      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 147      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 148      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 149      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 150      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 151      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 152      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 153      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 154      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 155
 156       if (wcorr4.gt.0.0d0) then
 157         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 158      &   'between contact pairs of peptide groups'
 159         write (iout,'(2(a,f5.3/))')
 160      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 161      &  'Range of quenching the correlation terms:',2*delt_corr
 162       else if (wcorr.gt.0.0d0) then
 163         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 164      &   'between contact pairs of peptide groups'
 165       endif
 166       write (iout,'(a,f8.3)')
 167      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 168       write (iout,'(a,f8.3)')
 169      &  'General scaling factor of SC-p interactions:',scalscp
 170       r0_corr=cutoff_corr-delt_corr
 171       do i=1,20
 172         aad(i,1)=scalscp*aad(i,1)
 173         aad(i,2)=scalscp*aad(i,2)
 174         bad(i,1)=scalscp*bad(i,1)
 175         bad(i,2)=scalscp*bad(i,2)
 176       enddo
 177
 178       call flush(iout)
 179       print *,'indpdb=',indpdb,' pdbref=',pdbref
 180
 181 C Read sequence if not taken from the pdb file.
 182       if (iscode.gt.0) then
 183         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 184       else
 185         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 186       endif
 187 C Convert sequence to numeric code
 188       do i=1,nres
 189         itype(i)=rescode(i,sequence(i),iscode)
 190       enddo
 191       print *,nres
 192       print '(20i4)',(itype(i),i=1,nres)
 193
 194       do i=1,nres
 195 #ifdef PROCOR
 196         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 197 #else
 198         if (itype(i).eq.ntyp1) then
 199 #endif
 200           itel(i)=0
 201 #ifdef PROCOR
 202         else if (iabs(itype(i+1)).ne.20) then
 203 #else
 204         else if (iabs(itype(i)).ne.20) then
 205 #endif
 206           itel(i)=1
 207         else
 208           itel(i)=2
 209         endif
 210       enddo
 211       write (iout,*) "ITEL"
 212       do i=1,nres-1
 213         write (iout,*) i,itype(i),itel(i)
 214       enddo
 215
 216       print *,'Call Read_Bridge.'
 217       call read_bridge
 218       nnt=1
 219       nct=nres
 220       print *,'NNT=',NNT,' NCT=',NCT
 221       if (itype(1).eq.ntyp1) nnt=2
 222       if (itype(nres).eq.ntyp1) nct=nct-1
 223       if (nstart.lt.nnt) nstart=nnt
 224       if (nend.gt.nct .or. nend.eq.0) nend=nct
 225       write (iout,*) "nstart",nstart," nend",nend
 226       nres0=nres
 227       nchain=1
 228       ichain(1,1)=nnt
 229       i=nnt
 230       nres_chain(nnt)=1
 231       do while (i.lt.nct)
 232         i=i+1
 233         if (itype(i).eq.ntyp1) then
 234           ichain(2,nchain)=i-1
 235           nchain=nchain+1
 236           do while (itype(i).eq.ntyp1)
 237             i=i+1
 238           enddo
 239           ichain(1,nchain)=i
 240         endif
 241         nres_chain(i)=nchain
 242       enddo
 243       ichain(2,nchain)=nct
 244       write(iout,*) "nchain",nchain
 245       do i=1,nchain
 246         write(iout,*) ichain(1,i),ichain(2,i)
 247       enddo
 248       call init_int_table
 249       call setup_var
 250       write (iout,*) "molread: REFSTR",refstr
 251       return
 252       end
 253 c-----------------------------------------------------------------------------
 254       logical function seq_comp(itypea,itypeb,length)
 255       implicit none
 256       integer length,itypea(length),itypeb(length)
 257       integer i
 258       do i=1,length
 259         if (itypea(i).ne.itypeb(i)) then
 260           seq_comp=.false.
 261           return
 262         endif
 263       enddo
 264       seq_comp=.true.
 265       return
 266       end
 267 c-----------------------------------------------------------------------------
 268       subroutine read_bridge
 269 C Read information about disulfide bridges.
 270       implicit none
 271       include 'DIMENSIONS'
 272       include 'COMMON.IOUNITS'
 273       include 'COMMON.GEO'
 274       include 'COMMON.VAR'
 275       include 'COMMON.INTERACT'
 276       include 'COMMON.LOCAL'
 277       include 'COMMON.NAMES'
 278       include 'COMMON.CHAIN'
 279       include 'COMMON.FFIELD'
 280       include 'COMMON.SBRIDGE'
 281       include 'COMMON.HEADER'
 282       include 'COMMON.CONTROL'
 283       include 'COMMON.TIME1'
 284 #ifdef MPL
 285       include 'COMMON.INFO'
 286 #endif
 287       integer i,j
 288 C Read bridging residues.
 289       read (inp,*) ns,(iss(i),i=1,ns)
 290       print *,'ns=',ns
 291 C Check whether the specified bridging residues are cystines.
 292       do i=1,ns
 293         if (itype(iss(i)).ne.1) then
 294           write (iout,'(2a,i3,a)')
 295      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 296      &   ' can form a disulfide bridge?!!!'
 297           write (*,'(2a,i3,a)')
 298      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 299      &   ' can form a disulfide bridge?!!!'
 300 #ifdef MPL
 301          call mp_stopall(error_msg)
 302 #else
 303          stop
 304 #endif
 305         endif
 306       enddo
 307 C Read preformed bridges.
 308       if (ns.gt.0) then
 309       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 310       if (nss.gt.0) then
 311         nhpb=nss
 312 C Check if the residues involved in bridges are in the specified list of
 313 C bridging residues.
 314         do i=1,nss
 315           do j=1,i-1
 316             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 317      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 318               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 319      &      ' contains residues present in other pairs.'
 320               write (*,'(a,i3,a)') 'Disulfide pair',i,
 321      &      ' contains residues present in other pairs.'
 322 #ifdef MPL
 323               call mp_stopall(error_msg)
 324 #else
 325               stop
 326 #endif
 327             endif
 328           enddo
 329           do j=1,ns
 330             if (ihpb(i).eq.iss(j)) goto 10
 331           enddo
 332           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 333    10     continue
 334           do j=1,ns
 335             if (jhpb(i).eq.iss(j)) goto 20
 336           enddo
 337           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 338    20     continue
 339           dhpb(i)=dbr
 340           forcon(i)=fbr
 341         enddo
 342         do i=1,nss
 343           ihpb(i)=ihpb(i)+nres
 344           jhpb(i)=jhpb(i)+nres
 345         enddo
 346       endif
 347       endif
 348       return
 349       end
 350 c----------------------------------------------------------------------------
 351       subroutine read_angles(kanal,*)
 352       implicit none
 353       include 'DIMENSIONS'
 354       include 'COMMON.GEO'
 355       include 'COMMON.VAR'
 356       include 'COMMON.CHAIN'
 357       include 'COMMON.IOUNITS'
 358       integer i,kanal
 359       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 360       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 361       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 362       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 363       do i=1,nres
 364         theta(i)=deg2rad*theta(i)
 365         phi(i)=deg2rad*phi(i)
 366         alph(i)=deg2rad*alph(i)
 367         omeg(i)=deg2rad*omeg(i)
 368       enddo
 369       return
 370    10 return1
 371       end
 372 c----------------------------------------------------------------------------
 373       subroutine reada(rekord,lancuch,wartosc,default)
 374       implicit none
 375       character*(*) rekord,lancuch
 376       double precision wartosc,default
 377       integer ilen,iread
 378       external ilen
 379       iread=index(rekord,lancuch)
 380       if (iread.eq.0) then
 381         wartosc=default
 382         return
 383       endif
 384       iread=iread+ilen(lancuch)+1
 385       read (rekord(iread:),*) wartosc
 386       return
 387       end
 388 c----------------------------------------------------------------------------
 389       subroutine multreada(rekord,lancuch,tablica,dim,default)
 390       implicit none
 391       integer dim,i
 392       double precision tablica(dim),default
 393       character*(*) rekord,lancuch
 394       integer ilen,iread
 395       external ilen
 396       do i=1,dim
 397         tablica(i)=default
 398       enddo
 399       iread=index(rekord,lancuch)
 400       if (iread.eq.0) return
 401       iread=iread+ilen(lancuch)+1
 402       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 403    10 return
 404       end
 405 c----------------------------------------------------------------------------
 406       subroutine readi(rekord,lancuch,wartosc,default)
 407       implicit none
 408       character*(*) rekord,lancuch
 409       integer wartosc,default
 410       integer ilen,iread
 411       external ilen
 412       iread=index(rekord,lancuch)
 413       if (iread.eq.0) then
 414         wartosc=default
 415         return
 416       endif
 417       iread=iread+ilen(lancuch)+1
 418       read (rekord(iread:),*) wartosc
 419       return
 420       end
 421 c----------------------------------------------------------------------------
 422       subroutine card_concat(card)
 423       include 'DIMENSIONS'
 424       include 'COMMON.IOUNITS'
 425       character*(*) card
 426       character*80 karta,ucase
 427       external ilen
 428       read (inp,'(a)') karta
 429       karta=ucase(karta)
 430       card=' '
 431       do while (karta(80:80).eq.'&')
 432         card=card(:ilen(card)+1)//karta(:79)
 433         read (inp,'(a)') karta
 434         karta=ucase(karta)
 435       enddo
 436       card=card(:ilen(card)+1)//karta
 437       return
 438       end
 439 c----------------------------------------------------------------------------
 440       subroutine openunits
 441       implicit none
 442       include 'DIMENSIONS'
 443 #ifdef MPI
 444       include "mpif.h"
 445       character*3 liczba
 446       include "COMMON.MPI"
 447 #endif
 448       include 'COMMON.IOUNITS'
 449       include 'COMMON.CONTROL'
 450       integer lenpre,lenpot,ilen
 451       external ilen
 452       character*16 cformat,cprint
 453       character*16 ucase
 454       integer lenint,lenout
 455       call getenv('INPUT',prefix)
 456       call getenv('OUTPUT',prefout)
 457       call getenv('INTIN',prefintin)
 458       call getenv('COORD',cformat)
 459       call getenv('PRINTCOOR',cprint)
 460       call getenv('SCRATCHDIR',scratchdir)
 461       from_bx=.true.
 462       from_cx=.false.
 463       if (index(ucase(cformat),'CX').gt.0) then
 464         from_cx=.true.
 465         from_bx=.false.
 466       endif
 467       from_cart=.true.
 468       lenpre=ilen(prefix)
 469       lenout=ilen(prefout)
 470       lenint=ilen(prefintin)
 471 C Get the names and open the input files
 472       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 473 #ifdef MPI
 474       write (liczba,'(bz,i3.3)') me
 475       outname=prefout(:lenout)//'_clust.out_'//liczba
 476 #else
 477       outname=prefout(:lenout)//'_clust.out'
 478 #endif
 479       if (from_bx) then
 480         intinname=prefintin(:lenint)//'.bx'
 481       else if (from_cx) then
 482         intinname=prefintin(:lenint)//'.cx'
 483       else
 484         intinname=prefintin(:lenint)//'.int'
 485       endif
 486       open(iout,file=outname,status='unknown')
 487       open(istat,file=prefix(:ilen(prefix))//'.stat',status='unknown')
 488 C Get parameter filenames and open the parameter files.
 489       call getenv('BONDPAR',bondname)
 490       open (ibond,file=bondname,status='old')
 491       call getenv('THETPAR',thetname)
 492       open (ithep,file=thetname,status='old')
 493       call getenv('ROTPAR',rotname)
 494       open (irotam,file=rotname,status='old')
 495       call getenv('TORPAR',torname)
 496       open (itorp,file=torname,status='old')
 497       call getenv('TORDPAR',tordname)
 498       open (itordp,file=tordname,status='old')
 499       call getenv('FOURIER',fouriername)
 500       open (ifourier,file=fouriername,status='old')
 501       call getenv('ELEPAR',elename)
 502       open (ielep,file=elename,status='old')
 503       call getenv('SIDEPAR',sidename)
 504       open (isidep,file=sidename,status='old')
 505       call getenv('SIDEP',sidepname)
 506       open (isidep1,file=sidepname,status="old")
 507       call getenv('SCCORPAR',sccorname)
 508       open (isccor,file=sccorname,status="old")
 509 #ifndef OLDSCP
 510 C
 511 C 8/9/01 In the newest version SCp interaction constants are read from a file
 512 C Use -DOLDSCP to use hard-coded constants instead.
 513 C
 514       call getenv('SCPPAR',scpname)
 515       open (iscpp,file=scpname,status='old')
 516 #endif
 517       return
 518       end