Add bpall file

Hori-Lab · Sep 12, 2021 · 90ab5e2 · 90ab5e2
1 parent d1cfb16
commit 90ab5e2
Show file tree

Hide file tree

Showing 9 changed files with 53 additions and 16 deletions.
diff --git a/input.toml b/input.toml
@@ -10,7 +10,7 @@ type = "DCD"
       fasta = "./beta-globin.fasta"
 
       [files.out]
-      types = ["bp", "bpe"]
+      types = ["bp", "bpe", "bpall"]
       prefix = "./test"
 
 [condition]

diff --git a/src/energy_bp.F90 b/src/energy_bp.F90
@@ -4,7 +4,7 @@ subroutine energy_bp(Ebp)
    use const_phys, only : ZERO_JUDGE
    use var_state, only : xyz, kT
    use var_potential
-   use var_io, only : flg_out_bp, flg_out_bpe, hdl_bp, hdl_bpe, KIND_OUT_BP, KIND_OUT_BPE
+   use var_io, only : flg_out_bp, flg_out_bpall, flg_out_bpe, hdl_bp, hdl_bpall, hdl_bpe, KIND_OUT_BP, KIND_OUT_BPE
 
    implicit none
 
@@ -64,15 +64,30 @@ subroutine energy_bp(Ebp)
          if (e_bp(ibp) < - nhb * kT) then
             imp = bp_mp(1, ibp)
             jmp = bp_mp(2, ibp)
-            !write(hdl_bp) int(imp,kind=KIND_OUT_BP), int(jmp,kind=KIND_OUT_BP)
             write(hdl_bp) int(imp,kind=KIND_OUT_BP), int(jmp,kind=KIND_OUT_BP), real(e_bp(ibp), kind=KIND_OUT_BPE)
          endif
       enddo
 
-      !write(hdl_bp) int(0,kind=KIND_OUT_BP)
       write(hdl_bp) int(0,kind=KIND_OUT_BP), int(0,kind=KIND_OUT_BP), real(0.0, kind=KIND_OUT_BPE)
    endif
 
+   if (flg_out_bpall) then
+
+      do ibp = 1, nbp
+
+         nhb = bp_mp(3, ibp)
+
+         if (e_bp(ibp) < -ZERO_JUDGE) then  ! To output all
+            imp = bp_mp(1, ibp)
+            jmp = bp_mp(2, ibp)
+            write(hdl_bpall) int(imp,kind=KIND_OUT_BP), int(jmp,kind=KIND_OUT_BP), real(e_bp(ibp), kind=KIND_OUT_BPE)
+         endif
+      enddo
+
+      write(hdl_bpall) int(0,kind=KIND_OUT_BP), int(0,kind=KIND_OUT_BP), real(0.0, kind=KIND_OUT_BPE)
+   endif
+
+
    if (flg_out_bpe) then
 
       do ibp = 1, nbp

diff --git a/src/force_bp.F90 b/src/force_bp.F90
@@ -4,7 +4,7 @@ subroutine force_bp()
    use const_phys, only : ZERO_JUDGE
    use var_state, only : xyz, kT, forces
    use var_potential
-   use var_io, only : flg_out_bp, flg_out_bpe, hdl_bp, hdl_bpe, KIND_OUT_BP
+   use var_io, only : flg_out_bp, flg_out_bpe, KIND_OUT_BP
 
    implicit none
 

diff --git a/src/main.F90 b/src/main.F90
@@ -6,13 +6,13 @@ program sis
    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_fasta_in
+   use var_io, only : flg_out_bp, flg_out_bpall, flg_out_bpe, hdl_out, hdl_bp, hdl_bpall, hdl_bpe, KIND_OUT_BP, KIND_OUT_BPE, &
+                      cfile_ff, cfile_dcd_in, cfile_prefix, cfile_out, cfile_fasta_in
 !$ use omp_lib
 
    implicit none
 
-   character(CHAR_FILE_PATH) cfile_inp
+   character(CHAR_FILE_PATH) :: cfile_inp, cfile_bp
 
    integer :: i, j, k, imp
    integer :: nthreads
@@ -63,8 +63,11 @@ program sis
       stop (2) 
    end if
 
-   write(*,*) 'Read force-field file: '//cfile_ff
-   call read_force_field(cfile_ff)
+   !! Load force field
+   write(*,*) 'Read force-field file: '//trim(cfile_ff)
+   call read_force_field(trim(cfile_ff))
+
+   !! Output files
    cfile_out = trim(cfile_prefix) // '.out'
 
    if (flg_out_bp) then
@@ -74,6 +77,13 @@ program sis
       write(hdl_bp) int(KIND_OUT_BPE,kind=4)
    endif
 
+   if (flg_out_bpall) then
+      cfile_bp = trim(cfile_prefix) // '.bpall'
+      open(hdl_bpall, file=cfile_bp, status='replace', action='write', form='unformatted',access='stream')
+      write(hdl_bpall) int(KIND_OUT_BP,kind=4)
+      write(hdl_bpall) int(KIND_OUT_BPE,kind=4)
+   endif
+
    if (flg_out_bpe) then
       cfile_bp = trim(cfile_prefix) // '.bpe'
       open(hdl_bpe, file=cfile_bp, status='replace', action='write', form='formatted')

diff --git a/src/read_fasta.F90 b/src/read_fasta.F90
@@ -138,5 +138,8 @@ subroutine read_fasta()
       endif
 
    enddo
+
+   close(hdl)
+   iopen_hdl = iopen_hdl - 1
 
 end subroutine read_fasta
diff --git a/src/read_force_field.F90 b/src/read_force_field.F90
@@ -14,9 +14,6 @@ subroutine read_force_field(cfilepath)
 
    integer :: hdl
 
-   iopen_hdl = iopen_hdl + 1
-   hdl = iopen_hdl
-
    bond_k  = INVALID_VALUE
    bond_r0 = INVALID_VALUE
    angl_k  = INVALID_VALUE
@@ -34,6 +31,9 @@ subroutine read_force_field(cfilepath)
    Ubp_dihd_phi2 = INVALID_VALUE
    Ubp_min_loop = -1
 
+   iopen_hdl = iopen_hdl + 1
+   hdl = iopen_hdl
+
    open(hdl, file=cfilepath, status='old', action='read', iostat=istat)
 
    do

diff --git a/src/read_input.F90 b/src/read_input.F90
@@ -6,7 +6,7 @@ subroutine read_input(cfilepath, stat)
    use const_phys
    use const_idx
    use var_io, only : iopen_hdl, &
-                      flg_out_bp, flg_out_bpe, &
+                      flg_out_bp, flg_out_bpe, flg_out_bpall, &
                       cfile_ff, cfile_dcd_in, &
                       cfile_prefix, cfile_pdb_ini, cfile_fasta_in
    use var_state, only : job, tempK
@@ -59,6 +59,9 @@ subroutine read_input(cfilepath, stat)
    else if (cline == 'CHECK_FORCE') then
       job = JOBT%CHECK_FORCE
 
+   else if (cline == 'MD') then
+      job = JOBT%MD
+
    else
       write(*,*) 'Unknown job type: '//trim(cline)
       stat = .False.
@@ -110,6 +113,8 @@ subroutine read_input(cfilepath, stat)
 
          if (cline == "bp") then
             flg_out_bp = .True.
+         else if (cline == "bpall") then
+            flg_out_bpall = .True.
          else if (cline == "bpe") then
             flg_out_bpe = .True.
          else

diff --git a/src/read_pdb.F90 b/src/read_pdb.F90
@@ -89,4 +89,7 @@ subroutine read_pdb(cfilepath, nmp, xyz)
    !nmp   = imp
    !nres  = ires
 
+   close(hdl)
+   iopen_hdl = iopen_hdl - 1
+
 end subroutine read_pdb
diff --git a/src/var_io.F90 b/src/var_io.F90
@@ -3,19 +3,20 @@ module var_io
    use const
 
    logical, save :: flg_out_bp = .False.
+   logical, save :: flg_out_bpall = .False.
    logical, save :: flg_out_bpe = .False.
    integer, parameter :: KIND_OUT_BP  = 2  ! Defines the format of bp output file.
    integer, parameter :: KIND_OUT_BPE = 4  ! Defines the format of bp output file.
 
    integer, parameter :: hdl_dcd = 20
    integer, parameter :: hdl_out = 21
    integer, parameter :: hdl_bp  = 22
-   integer, parameter :: hdl_bpe = 23
+   integer, parameter :: hdl_bpall = 23
+   integer, parameter :: hdl_bpe = 24
 
    integer, save :: iopen_hdl = 30
 
    character(CHAR_FILE_PATH), save :: cfile_out
-   character(CHAR_FILE_PATH), save :: cfile_bp
    character(len=:), allocatable, save :: cfile_ff, cfile_prefix, cfile_dcd_in, cfile_pdb_ini, cfile_fasta_in
 
 end module var_io