Adding MARTINI

[unres4.git] / source / unres / io_base.F90

diff --git a/source/unres/io_base.F90 b/source/unres/io_base.F90
index e82c308..198fb1f 100644 (file)
--- a/source/unres/io_base.F90
+++ b/source/unres/io_base.F90
@@ -958,8 +958,9 @@
 !-----------------------------------------------------------------------------
       subroutine pdbout(etot,tytul,iunit)
 
-      use geometry_data, only: c,nres
+      use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize
       use energy_data
+!      use energy, only: to_box
 !     use control
       use compare_data
       use MD_data
@@ -975,7 +976,10 @@
 !      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(58) :: chainid=  (/'A','B','C','D','E','F','G','H','I','J',&
+         'K','L','M','O','Q','P','R','S','T','U','W','X','Y','Z','a','b','c','d','e','f',&
+            'g','h','i','j','k','l','m','n','o','p','r','s','t','u','w','x','y','z',&
+            '0','1','2','3','4','5','6','7','8','9'/)
 #ifdef XDRFPDB
       integer,dimension(maxres) :: ica !(maxres)
 #else
@@ -984,13 +988,14 @@
       integer :: iti1
 !el  local variables
       integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
-      real(kind=8) :: etot
+      real(kind=8) :: etot,xi,yi,zi
       integer :: nres2
       nres2=2*nres
 #ifdef XDRFPDB
       if(.not.allocated(molnum)) allocate(molnum(maxres))
-      molnum(:)=1
+!      molnum(:)=mnumi(:)
       if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2)
+      if(.not.allocated(istype)) allocate(istype(maxres))
 #else
       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
 #endif
@@ -1091,6 +1096,7 @@
         endif
         if (iti.eq.ntyp1_molec(molnum(i))) then
           ichain=ichain+1
+          if (ichain.gt.58) ichain=1
 !          ires=0
           write (iunit,'(a)') 'TER'
         else
@@ -1105,7 +1111,19 @@
            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
+        elseif (molnum(i).eq.4) then
+        xi=c(1,i)
+        yi=c(2,i)
+        zi=c(3,i)
+        xi=dmod(xi,boxxsize)
+        if (xi.lt.0.0d0) xi=xi+boxxsize
+        yi=dmod(yi,boxysize)
+        if (yi.lt.0.0d0) yi=yi+boxysize
+        zi=dmod(zi,boxzsize)
+        if (zi.lt.0.0d0) zi=zi+boxzsize
+        write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
+           ires,xi,yi,zi,vtot(i)
+        else 
         write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
            ires,(c(j,i),j=1,3),vtot(i)
         endif
@@ -1521,6 +1539,7 @@
       else if (molecule.eq.2) then
       do i=1,ntyp1_molec(molecule)
          print *,nam(1:1),restyp(i,molecule)(1:1)
+         print *,nam(2:2),"read",i
         if (nam(2:2).eq.restyp(i,molecule)(1:1)) then
           rescode=i
           return
@@ -1561,6 +1580,29 @@
       endif
       return
       end function sugarcode
+      integer function rescode_lip(res,ires)
+      character(len=3):: res
+      integer ires
+       rescode_lip=0
+       if  (res.eq.'NC3')  rescode_lip =4
+       if  (res.eq.'NH3')  rescode_lip =2
+       if  (res.eq.'PO4')  rescode_lip =3
+       if  (res.eq.'GL1')  rescode_lip =12
+       if  (res.eq.'GL2')  rescode_lip =12
+       if  (res.eq.'C1A')  rescode_lip =18
+       if  (res.eq.'C2A')  rescode_lip =18
+       if  (res.eq.'C3A')  rescode_lip =18
+       if  (res.eq.'C4A')  rescode_lip =18
+       if  (res.eq.'C1B')  rescode_lip =18
+       if  (res.eq.'C2B')  rescode_lip =18
+       if  (res.eq.'C3B')  rescode_lip =18
+       if  (res.eq.'C4B')  rescode_lip =18
+       if  (res.eq.'D2A')  rescode_lip =16
+       if  (res.eq.'D2B')  rescode_lip =16
+
+        if (rescode_lip.eq.0) write(iout,*) "UNKNOWN RESIDUE",ires,res
+      return
+      end function
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
       end module io_base