uses only first alternate conformation, and prints warning accordingly

m3g · Jun 30, 2018 · 5bb34f7 · 5bb34f7
1 parent 4d4dd7e
commit 5bb34f7
Show file tree

Hide file tree

Showing 3 changed files with 93 additions and 6 deletions.
diff --git a/src/Makefile b/src/Makefile
@@ -66,7 +66,6 @@ $(BIN)/topolink: Makefile $(MODULES) \
                  setskip.o \
                  solventaccess.o \
                  setsigma.o \
-                 isprotein.o \
                  test_grad.o \
                  title.o
 	$(FORTRAN) -o $(BIN)/topolink $(MODULES) \
@@ -91,7 +90,6 @@ $(BIN)/topolink: Makefile $(MODULES) \
                  setskip.o \
                  solventaccess.o \
                  setsigma.o \
-                 isprotein.o \
                  test_grad.o \
                  title.o \
                  $(FLAGS)
@@ -227,9 +225,6 @@ solventaccess.o : solventaccess.f90
 setsigma.o : setsigma.f90
 	$(FORTRAN) -c setsigma.f90 $(FLAGS)
 
-isprotein.o : isprotein.f90
-	$(FORTRAN) -c isprotein.f90 $(FLAGS)
-
 test_grad.o : test_grad.f90
 	$(FORTRAN) -c test_grad.f90 $(FLAGS)
 

diff --git a/src/topolink.f90 b/src/topolink.f90
@@ -36,7 +36,7 @@ program topolink
   character(len=max_string_length), allocatable :: logline(:)
   character(len=20) :: floatout, intout, intout2
   character(len=13) :: statuschar
-  logical :: error, r1, r2, inexp, warning, interchain, isprotein
+  logical :: error, r1, r2, inexp, warning, interchain
 
   external :: computef, computeg
 
@@ -528,6 +528,14 @@ program topolink
     readatom=read_atom(record,error)
     if ( error ) cycle
     if ( .not. isprotein(readatom) ) cycle
+    if ( readatom%residue%conformation /= " " .and. &
+         readatom%residue%conformation /= "A" ) then
+      write(*,"(3(a,tr1),i8,a)") '  WARNING: Multiple conformations of residue',&
+                                 trim(adjustl(readatom%residue%name)),&
+                                 trim(adjustl(readatom%residue%chain)),&
+                                 readatom%residue%index,&
+                                 ' found. Using only A.'
+    end if
     if ( readatoms == 1 ) then
       if ( ishydrogen(readatom) ) cycle
     else if ( readatoms == 2 ) then

diff --git a/src/topolink_data.f90 b/src/topolink_data.f90
@@ -13,6 +13,7 @@ module topolink_data
 
     integer :: index, firstatom, lastatom
     character(len=4) :: name, chain = "    "
+    character(len=1) :: conformation
     logical :: accessible
 
   end type pdbresidue
@@ -177,6 +178,7 @@ function read_atom(record,error)
          read(record(17:21),*,iostat=ioerr) read_atom%residue%name
          read_atom%residue%name = trim(adjustl(read_atom%residue%name))
          if ( ioerr /= 0 ) error = .true.
+         call alternate_conformation(read_atom%residue)
 
          if ( record(22:22) /= " " ) then
            read(record(22:22),*,iostat=ioerr) read_atom%residue%chain
@@ -428,5 +430,87 @@ function print_pdbhetatm(atom)
 
      end function print_pdbhetatm
 
+     !
+     ! function isprotein: check if an atom is a protein atom
+     !
+     logical function isprotein(atom)
+
+       implicit none
+       type(pdbatom) :: atom
+
+       isprotein = .false.
+       select case ( atom%residue%name ) 
+         case ( "ALA" ) ; isprotein = .true. ; return
+         case ( "ARG" ) ; isprotein = .true. ; return
+         case ( "ASN" ) ; isprotein = .true. ; return
+         case ( "ASP" ) ; isprotein = .true. ; return
+         case ( "ASX" ) ; isprotein = .true. ; return
+         case ( "CYS" ) ; isprotein = .true. ; return
+         case ( "GLU" ) ; isprotein = .true. ; return
+         case ( "GLN" ) ; isprotein = .true. ; return
+         case ( "GLX" ) ; isprotein = .true. ; return
+         case ( "GLY" ) ; isprotein = .true. ; return
+         case ( "HIS" ) ; isprotein = .true. ; return
+         case ( "HSE" ) ; isprotein = .true. ; return
+         case ( "HSD" ) ; isprotein = .true. ; return
+         case ( "ILE" ) ; isprotein = .true. ; return
+         case ( "LEU" ) ; isprotein = .true. ; return
+         case ( "LYS" ) ; isprotein = .true. ; return
+         case ( "MET" ) ; isprotein = .true. ; return
+         case ( "PHE" ) ; isprotein = .true. ; return
+         case ( "PRO" ) ; isprotein = .true. ; return
+         case ( "SER" ) ; isprotein = .true. ; return
+         case ( "THR" ) ; isprotein = .true. ; return
+         case ( "TRP" ) ; isprotein = .true. ; return
+         case ( "TYR" ) ; isprotein = .true. ; return
+         case ( "VAL" ) ; isprotein = .true. ; return
+         case default ; return
+       end select
+
+     end function isprotein
+
+     !
+     ! subroutine alternate_conformation: check if an atom is from an alternate
+     !                                    conformation, and pick only the "A" conformation
+     !
+     !
+     subroutine alternate_conformation(residue)
+
+       implicit none
+       type(pdbresidue) :: residue
+
+       select case ( residue%name(2:4) ) 
+         case ( "ALA" , &
+                "ARG" , & 
+                "ASN" , & 
+                "ASP" , & 
+                "ASX" , & 
+                "CYS" , & 
+                "GLU" , & 
+                "GLN" , & 
+                "GLX" , & 
+                "GLY" , & 
+                "HIS" , & 
+                "HSE" , & 
+                "HSD" , & 
+                "ILE" , & 
+                "LEU" , & 
+                "LYS" , & 
+                "MET" , & 
+                "PHE" , & 
+                "PRO" , & 
+                "SER" , & 
+                "THR" , & 
+                "TRP" , & 
+                "TYR" , & 
+                "VAL" ) 
+           residue%conformation = residue%name(1:1)
+           residue%name = residue%name(2:4)
+         case default 
+           residue%conformation = " "
+       end select
+
+     end subroutine alternate_conformation
+
 end module topolink_data