X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_base.f90;h=d07c1f00612d5923ff3c03e4e34acdbf1ccbeb41;hb=10689ab7d813dfbdbb0c6e631d90234d78ea306a;hp=dd8d9ed7228eeb23ecc189667b7a0c5a551a56ba;hpb=c50b0dd5a791f60040b970bcaf655751d286dce1;p=unres4.git diff --git a/source/unres/io_base.f90 b/source/unres/io_base.f90 index dd8d9ed..d07c1f0 100644 --- a/source/unres/io_base.f90 +++ b/source/unres/io_base.f90 @@ -5,7 +5,7 @@ implicit none !----------------------------------------------------------------------------- ! Max. number of AA residues - integer,parameter :: maxres=6000!1200 + integer,parameter :: maxres=6500!1200 ! Appr. max. number of interaction sites integer,parameter :: maxres2=2*maxres ! parameter (maxres6=6*maxres) @@ -25,6 +25,7 @@ !----------------------------------------------------------------------------- ! readrtns_CSA.F !----------------------------------------------------------------------------- +#ifndef XDRFPDB subroutine read_bridge ! Read information about disulfide bridges. use geometry_data, only: nres @@ -56,7 +57,8 @@ print *,"ENTER READ" ! Read bridging residues. read (inp,*) ns - print *,"ns",ns + write(iout,*) "ns",ns + call flush(iout) if (ns.gt.0) then allocate(iss(ns)) read (inp,*) (iss(i),i=1,ns) @@ -275,8 +277,9 @@ subroutine read_angles(kanal,*) use geometry_data - ! use energy - ! use control + +! use energy +! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.GEO' @@ -452,7 +455,7 @@ !----------------------------------------------------------------------------- subroutine read_dist_constr use MPI_data - ! use control +! use control use geometry, only: dist use geometry_data use control_data @@ -485,6 +488,7 @@ if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim)) if(.not.allocated(forcon)) allocate(forcon(maxdim)) if(.not.allocated(fordepth)) allocate(fordepth(maxdim)) + if(.not.allocated(ibecarb)) allocate(ibecarb(maxdim)) if ((genconstr.gt.0).and.(constr_dist.eq.11)) then call gen_dist_constr2 go to 1712 @@ -593,7 +597,9 @@ if (constr_dist.eq.11) then read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) - fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) +!EL fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) + fordepth(nhpb+1)=fordepth(nhpb+1)**(0.25d0) + forcon(nhpb+1)=forcon(nhpb+1)**(0.25d0) else !C print *,"in else" read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & @@ -625,7 +631,7 @@ !----------------------------------------------------------------------------- subroutine gen_dist_constr2 use MPI_data - ! use control +! use control use geometry, only: dist use geometry_data use control_data @@ -772,6 +778,7 @@ 10 find_group=.false. return end function find_group +#endif !----------------------------------------------------------------------------- logical function iblnk(charc) character(len=1) :: charc @@ -781,6 +788,7 @@ return end function iblnk !----------------------------------------------------------------------------- +#ifndef XDRFPDB integer function ilen(string) character*(*) :: string !EL logical :: iblnk @@ -918,6 +926,7 @@ numm(1:1)=huj(inum2+1:inum2+1) return end subroutine numstr +#endif !----------------------------------------------------------------------------- function ucase(string) integer :: i, k, idiff @@ -951,7 +960,7 @@ use geometry_data, only: c,nres use energy_data - ! use control +! use control use compare_data use MD_data ! implicit real*8 (a-h,o-z) @@ -966,19 +975,29 @@ ! include 'COMMON.MD' !el character(len=50) :: tytul character*(*) :: tytul - character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/) + character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/) +#ifdef XDRFPDB + integer,dimension(maxres) :: ica !(maxres) +#else integer,dimension(nres) :: ica !(maxres) - integer iti1 +#endif + integer :: iti1 !el local variables integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit real(kind=8) :: etot integer :: nres2 nres2=2*nres - +#ifdef XDRFPDB + if(.not.allocated(molnum)) allocate(molnum(maxres)) + molnum(:)=1 + if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2) +#else if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2) +#endif write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot !model write (iunit,'(a5,i6)') 'MODEL',1 +#ifndef XDRFPDB if (nhfrag.gt.0) then do j=1,nhfrag iti=itype(hfrag(1,j),1) @@ -1040,7 +1059,7 @@ enddo endif - +#endif if (nss.gt.0) then do i=1,nss if (dyn_ss) then @@ -1070,7 +1089,7 @@ if (i.lt.nnt) then if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle endif - if (iti.eq.ntyp1) then + if (iti.eq.ntyp1_molec(molnum(i))) then ichain=ichain+1 ! ires=0 write (iunit,'(a)') 'TER' @@ -1078,18 +1097,30 @@ ires=ires+1 iatom=iatom+1 ica(i)=iatom - if (molnum(i).ne.5) then + if (molnum(i).eq.1) then + write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,(c(j,i),j=1,3),vtot(i) + elseif(molnum(i).eq.2) then + if (istype(i).eq.0) istype(i)=1 + write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), & + chainid(ichain),ires,(c(j,i),j=1,3),vtot(i) else write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,(c(j,i),j=1,3),vtot(i) endif if ((iti.ne.10).and.(molnum(i).ne.5)) then iatom=iatom+1 - write (iunit,20) iatom,restyp(iti,1),chainid(ichain),& + if (molnum(i).eq.1) then + write (iunit,20) iatom,restyp(iti,1),chainid(ichain),& ires,(c(j,nres+i),j=1,3),& vtot(i+nres) + else if (molnum(i).eq.2) then + if (istype(i).eq.0) istype(i)=1 + write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), & + chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres) + endif + endif endif enddo @@ -1118,18 +1149,22 @@ 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3) 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3) + 40 FORMAT ("ATOM",I7," C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3) + 50 FORMAT ("ATOM",I7," C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3) + 30 FORMAT ('CONECT',8I5) 60 FORMAT ('HETATM',I5,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3) return end subroutine pdbout +#ifndef XDRFPDB !----------------------------------------------------------------------------- subroutine MOL2out(etot,tytul) ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2 ! format. use geometry_data, only: c use energy_data - ! use control +! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' @@ -1193,7 +1228,7 @@ use geometry_data use energy_data, only: itype - ! use control +! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' @@ -1227,7 +1262,7 @@ #endif use geometry_data use energy_data - ! use control +! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' @@ -1296,7 +1331,7 @@ #endif use geometry_data, only: c - ! use geometry +! use geometry use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -1448,6 +1483,83 @@ nextp=.true. return end function nextp +#endif +!----------------------------------------------------------------------------- +! rescode.f +!----------------------------------------------------------------------------- + integer function rescode(iseq,nam,itype,molecule) + +! use io_base, only: ucase +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.NAMES' +! include 'COMMON.IOUNITS' + character(len=3) :: nam !,ucase + integer :: iseq,itype,i + integer :: molecule + print *,molecule,nam + if (molecule.eq.1) then + if (itype.eq.0) then + + do i=-ntyp1_molec(molecule),ntyp1_molec(molecule) + if (ucase(nam).eq.restyp(i,molecule)) then + rescode=i + return + endif + enddo + + else + + do i=-ntyp1_molec(molecule),ntyp1_molec(molecule) + if (nam(1:1).eq.onelet(i)) then + rescode=i + return + endif + enddo + + endif + else if (molecule.eq.2) then + do i=1,ntyp1_molec(molecule) + print *,nam(1:1),restyp(i,molecule)(1:1) + if (nam(2:2).eq.restyp(i,molecule)(1:1)) then + rescode=i + return + endif + enddo + else if (molecule.eq.3) then + write(iout,*) "SUGAR not yet implemented" + stop + else if (molecule.eq.4) then + write(iout,*) "Explicit LIPID not yet implemented" + stop + else if (molecule.eq.5) then + do i=1,ntyp1_molec(molecule) + print *,i,restyp(i,molecule)(1:2) + if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then + rescode=i + return + endif + enddo + else + write(iout,*) "molecule not defined" + endif + write (iout,10) iseq,nam + stop + 10 format ('**** Error - residue',i4,' has an unresolved name ',a3) + end function rescode + integer function sugarcode(sugar,ires) + character sugar + integer ires + if (sugar.eq.'D') then + sugarcode=1 + else if (sugar.eq.' ') then + sugarcode=2 + else + write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar + stop + endif + return + end function sugarcode !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module io_base