Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / unres / io_base.f90
index dd8d9ed..d07c1f0 100644 (file)
@@ -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)
       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'
 !-----------------------------------------------------------------------------
       subroutine read_dist_constr
       use MPI_data
-   !  use control
+!     use control
       use geometry, only: dist
       use geometry_data
       use control_data
       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
         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), &
 !-----------------------------------------------------------------------------
       subroutine gen_dist_constr2
       use MPI_data
-   !  use control
+!     use control
       use geometry, only: dist
       use geometry_data
       use control_data
    10 find_group=.false.
       return
       end function find_group
+#endif
 !-----------------------------------------------------------------------------
       logical function iblnk(charc)
       character(len=1) :: charc
       return
       end function iblnk
 !-----------------------------------------------------------------------------
+#ifndef XDRFPDB
       integer function ilen(string)
       character*(*) ::  string
 !EL      logical :: iblnk
       numm(1:1)=huj(inum2+1:inum2+1)
       return
       end subroutine numstr
+#endif
 !-----------------------------------------------------------------------------
       function ucase(string)
       integer :: i, k, idiff
 
       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)
 !      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)
          
        enddo
       endif 
-
+#endif
       if (nss.gt.0) then
         do i=1,nss
          if (dyn_ss) then
         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'
         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
   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'
 
       use geometry_data
       use energy_data, only: itype
-   !  use control
+!      use control
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.IOUNITS'
 #endif
       use geometry_data
       use energy_data
-   !  use control
+!      use control
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.IOUNITS'
 #endif
 
       use geometry_data, only: c
-   !  use geometry
+!      use geometry
       use energy_data
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
       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