Add FASTA format reader

input.toml

@@ -5,14 +5,10 @@ type = "DCD"
 
 [files]
       [files.in]
-      ff = "./rna_cg1.ff"
-      dcd = "md.dcd"
+      ff = "./rna_cg1_para2.ff"
+      dcd = "./HBB2.dcd"
+      fasta = "./beta-globin.fasta"
 
       [files.out]
       types = ["bp", "bpe"]
       prefix = "./test"
-
-
-[repeat]
-n_repeat = 47
-n_chain = 64

input_repeat.toml

@@ -0,0 +1,18 @@
+title = "sis input"
+
+[job]
+type = "DCD"
+
+[files]
+      [files.in]
+      ff = "./rna_cg1.ff"
+      dcd = "md.dcd"
+
+      [files.out]
+      types = ["bp", "bpe"]
+      prefix = "./test"
+
+
+[repeat]
+n_repeat = 47
+n_chain = 64

src/CMakeLists.txt

@@ -3,7 +3,7 @@ set(SIS_EXE sis)
 add_executable(${SIS_EXE}
     const.F90 const_idx.F90 const_phys.F90
     var_top.F90 var_state.F90 var_potential.F90 var_io.F90
-    dcd.F90 read_force_field.F90 read_input.F90 read_pdb.F90
+    dcd.F90 read_force_field.F90 read_input.F90 read_pdb.F90 read_fasta.F90
     list_local.F90 list_bp.F90
     energy.F90
     energy_bond.F90 energy_angl.F90 energy_bp.F90 read_sisinfo.F90

src/const_idx.F90

@@ -1,36 +1,81 @@
 module const_idx
 
-  use,intrinsic :: ISO_FORTRAN_ENV, only: REAL64, INT64
-  use :: const
-
-  implicit none
-
-  type energy_types
-     integer :: TOTAL  ! 0
-     integer :: BOND   ! 1
-     integer :: ANGL   ! 2
-     integer :: BP     ! 3
-     integer :: ELE    ! 4
-     integer :: MAX
-  endtype energy_types
-  type(energy_types), parameter :: ENE = energy_types(0,1,2,3,4,4)
-
-  type seq_types
-     integer :: A   ! 1
-     integer :: U   ! 2
-     integer :: G   ! 3
-     integer :: C   ! 4
-     integer :: MAX
-  endtype seq_types
-  type(seq_types), parameter :: SEQT = seq_types(1,2,3,4,4)
-
-  type job_types
-     integer :: DEBUG       ! 0
-     integer :: CHECK_FORCE ! 1
-     integer :: MD          ! 2
-     integer :: DCD         ! 3
-     integer :: MAX
-  endtype job_types
-  type(job_types), parameter :: JOBT = job_types(0,1,2,3,3)
+   use,intrinsic :: ISO_FORTRAN_ENV, only: REAL64, INT64
+   use :: const
+
+   implicit none
+
+   type energy_types
+      integer :: TOTAL  ! 0
+      integer :: BOND   ! 1
+      integer :: ANGL   ! 2
+      integer :: BP     ! 3
+      integer :: ELE    ! 4
+      integer :: MAX
+   endtype energy_types
+   type(energy_types), parameter :: ENE = energy_types(0,1,2,3,4,4)
+
+   type seq_types
+      integer :: UNDEF  ! 0
+      integer :: A   ! 1
+      integer :: U   ! 2
+      integer :: G   ! 3
+      integer :: C   ! 4
+      integer :: MAX
+   endtype seq_types
+   type(seq_types), parameter :: SEQT = seq_types(0,1,2,3,4,4)
+
+   type job_types
+      integer :: DEBUG       ! 0
+      integer :: CHECK_FORCE ! 1
+      integer :: MD          ! 2
+      integer :: DCD         ! 3
+      integer :: MAX
+   endtype job_types
+   type(job_types), parameter :: JOBT = job_types(0,1,2,3,3)
+
+contains
+   function seqt2char(i) result(c)
+      integer, intent(in) :: i
+      character(1) :: c
+
+      select case (i)
+      case (SEQT%A)
+         c = 'A'
+      case (SEQT%U)
+         c = 'U'
+      case (SEQT%G)
+         c = 'G'
+      case (SEQT%C)
+         c = 'C'
+      case (SEQT%UNDEF)
+         c = '?'
+      case default
+         c = '?'
+      endselect
+
+   endfunction seqt2char
+
+   function char2seqt(c) result (i)
+      character(1), intent(in) :: c
+      integer :: i
+
+      if (c == 'A' .or. c == 'a') then
+         i = SEQT%A
+
+      else if (c == 'U' .or. c == 'u') then
+         i = SEQT%U
+
+      else if (c == 'G' .or. c == 'g') then
+         i = SEQT%G
+
+      else if (c == 'C' .or. c == 'c') then
+         i = SEQT%C
+
+      else
+         i = SEQT%UNDEF
+      endif
+
+   endfunction char2seqt
 
 end module const_idx

src/main.F90

@@ -3,18 +3,18 @@ program sis
    use, intrinsic :: iso_fortran_env, Only : iostat_end
    use const
    use const_phys, only : BOLTZ_KCAL_MOL
-   use const_idx, only : ENE, SEQT, JOBT
+   use const_idx, only : ENE, SEQT, JOBT, seqt2char
    use var_top, only : nmp, nchains, nmp_chain, seq, imp_chain, pbc_box, pbc_box_half, flg_pbc, ichain_mp, nrepeat
    use var_state, only : xyz, tempK, kT, job
    use var_io, only : flg_out_bp, flg_out_bpe, hdl_out, hdl_bp, hdl_bpe, KIND_OUT_BP, KIND_OUT_BPE, &
-                      cfile_ff, cfile_dcd_in, cfile_prefix, cfile_out, cfile_bp
+                      cfile_ff, cfile_dcd_in, cfile_prefix, cfile_out, cfile_bp, cfile_fasta_in
 !$ use omp_lib
 
    implicit none
 
    character(CHAR_FILE_PATH) cfile_inp
 
-   integer :: i, j, imp
+   integer :: i, j, k, imp
    integer :: nthreads
 
    character(500) :: cline
@@ -80,31 +80,63 @@ program sis
 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-   allocate(nmp_chain(nchains))
-   nmp_chain(:) = 3 * nrepeat
-   nmp = sum(nmp_chain)
-   write(*, '(a,i8)') '#Nchain: ', nchains
-   write(*, '(a,i5)') '#Nrepeat: ', nrepeat
-   write(*, '(a,i10)') '#Nnt: ', nmp
-   allocate(seq(3*nrepeat, nchains))
-   allocate(imp_chain(3*nrepeat, nchains))
-   allocate(ichain_mp(nmp))
-   imp = 0
+   if (allocated(cfile_fasta_in)) then
+
+      write(*,*) 'reading fasta'
+      call read_fasta()
+
+   else if (nrepeat > 0) then
+
+      allocate(nmp_chain(nchains))
+      nmp_chain(:) = 3 * nrepeat
+      nmp = sum(nmp_chain)
+      write(*, '(a,i5)') '#Nrepeat: ', nrepeat
+      allocate(seq(3*nrepeat, nchains))
+      allocate(imp_chain(3*nrepeat, nchains))
+      allocate(ichain_mp(nmp))
+      imp = 0
+      do i = 1, nchains
+         do j = 1, nrepeat
+            seq(3*(j-1)+1, i) = SEQT%C
+            seq(3*(j-1)+2, i) = SEQT%A
+            seq(3*(j-1)+3, i) = SEQT%G
+            imp_chain(3*(j-1)+1, i) = imp+1
+            imp_chain(3*(j-1)+2, i) = imp+2
+            imp_chain(3*(j-1)+3, i) = imp+3
+            ichain_mp(imp+1) = i
+            ichain_mp(imp+2) = i
+            ichain_mp(imp+3) = i
+            imp = imp + 3
+         enddo
+         !write(*,'(a,141(i1))') '# ', seq(:,i)
+         !write(*,*) '# ', imp_chain(1,i), imp_chain((nrepeat-1)*3+3, i)
+      enddo
+
+   else
+      write(6,*) 'Error: either FASTA or [repeat] is required.'
+      stop (2)
+
+   endif
+
+   write(6, '(a)') '########### System'
+   write(6, '(a,i8)') '# Nchain: ', nchains
+   write(6, '(a)') '#'
    do i = 1, nchains
-      do j = 1, nrepeat
-         seq(3*(j-1)+1, i) = SEQT%C
-         seq(3*(j-1)+2, i) = SEQT%A
-         seq(3*(j-1)+3, i) = SEQT%G
-         imp_chain(3*(j-1)+1, i) = imp+1
-         imp_chain(3*(j-1)+2, i) = imp+2
-         imp_chain(3*(j-1)+3, i) = imp+3
-         ichain_mp(imp+1) = i
-         ichain_mp(imp+2) = i
-         ichain_mp(imp+3) = i
-         imp = imp + 3
+      write(6, '(a, i4)') '# Chain ', i
+      write(6, '(a, i10)') '# Nnt: ', nmp_chain(i)
+      k = 0
+      do j = 1, nmp_chain(i)
+         write(6, '(a)', advance='no') seqt2char(seq(j,i))
+         k = k + 1
+         if (mod(k,100) == 0) then
+            write(6, *) ''
+            k = 0
+         endif
       enddo
-      !write(*,'(a,141(i1))') '# ', seq(:,i)
-      !write(*,*) '# ', imp_chain(1,i), imp_chain((nrepeat-1)*3+3, i)
+      if (k /= 0) then
+         write(6, *) ''
+      endif
+      write(6, '(a)') '#'
    enddo
 
    tempK = 273.15 + 22.0