Add electrostatic (DH) potential

Hori-Lab · May 4, 2022 · 6ec3fed · 6ec3fed
1 parent ee5d794
commit 6ec3fed
Show file tree

Hide file tree

Showing 20 changed files with 542 additions and 64 deletions.
diff --git a/input_all.toml b/input_all.toml
@@ -49,6 +49,7 @@ x = 809.972
 y = 809.972
 z = 809.972
 
+
 [MD]
 integrator = "GJF-2GJ"
 dt = 0.05
@@ -60,9 +61,29 @@ neighbor_list_margin = 10.0
 viscosity_Pas = 0.00001
    # (optional) default is 0.00001
 
+
+[Electrostatic]
+ionic_strength = 0.1
+   # Ionic strength of the monovalent-ion solution
+
+cutoff_type = 1
+cutoff = 50.0
+   # cutoff_type: How to specify the cutoff distance for electrostatic interactions.
+   #             = 1: cutoff will be specified as distance in Angstrom. (default)
+   #             = 2: cutoff will be a factor to be multiplied by Debye length.
+   # cutoff: Either diestance (1) or factor (2), depending on the choice of cutoff_type.
+
+no_charge = [1, 27]
+   # Particles having no charges. This 
+
+length_per_charge = 4.38178046
+   # Paremeter in ion-condensation theory.
+
+
 [Progress]
 step = 1000
 
+
 [variable_box]
 step = 1000
 change_x = -1.0

diff --git a/rna_cg1.ff b/rna_cg1.ff
@@ -12,7 +12,9 @@ title = 'rna_cg1'
     [potential.basepair]
         min_loop = 4
         cutoff = 18.0
-        U0 = -1.66666666667
+        U0_GC = -5.00
+        U0_AU = -3.33
+        U0_GU = -3.33
 
         bond_k = 3.0
         bond_r = 13.8

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -9,14 +9,15 @@ add_executable(${SIS_EXE}
     read_force_field.F90 read_input.F90 read_pdb.F90 read_xyz.F90 read_fasta.F90
     read_anneal.F90
     write_rst.F90 read_rst.F90
-    list_local.F90 list_bp.F90 list_exv.F90
+    list_local.F90 list_bp.F90 list_exv.F90 list_ele.F90
+    set_ele.F90
     energy.F90
     energy_kinetic.F90
     energy_bond.F90 energy_angl.F90 energy_bp.F90 read_sisinfo.F90
-    energy_wca.F90
+    energy_wca.F90 energy_ele_DH.F90
     force.F90
     force_bond.F90 force_angl.F90 force_bp.F90
-    force_wca.F90
+    force_wca.F90 force_ele_DH.F90
     neighbor_list.F90
     job_md.F90 job_dcd.F90 job_check_force.F90
     main.F90

diff --git a/src/const_idx.F90 b/src/const_idx.F90
@@ -11,9 +11,10 @@ module const_idx
       integer :: ANGL   ! 2
       integer :: BP     ! 3
       integer :: EXV    ! 4
+      integer :: ELE    ! 5
       integer :: MAX
    endtype energy_types
-   type(energy_types), parameter :: ENE = energy_types(0,1,2,3,4,4)
+   type(energy_types), parameter :: ENE = energy_types(0,1,2,3,4,5,5)
 
    type seq_types
       integer :: UNDEF  ! 0

diff --git a/src/const_phys.F90 b/src/const_phys.F90
@@ -3,7 +3,7 @@ module const_phys
    use const
 
    real(PREC), parameter :: PI = 3.141592653589793238462643383279
-   real(PREC), parameter :: EPS0 = 8.8541878128             ! Electric constant [F/m]
+   real(PREC), parameter :: EPS0 = 8.8541878128e-12_PREC    ! Electric constant [F/m]
    real(PREC), parameter :: ELE  = 1.602176634e-19_PREC     ! Elementary charge [C]
    real(PREC), parameter :: BOLTZ_J = 1.380649e-23_PREC     ! Boltzmann constant [J/K]
    real(PREC), parameter :: N_AVO   = 6.02214076e23_PREC    ! Avogadro constant [/mol]

diff --git a/src/energy.F90 b/src/energy.F90
@@ -3,6 +3,7 @@ subroutine energy()
    use const
    use const_idx, only : ENE
    use var_state, only : energies
+   use var_potential, only : flg_ele
 
    implicit none
 
@@ -13,6 +14,8 @@ subroutine energy()
    call energy_bp(energies(ENE%BP))
    call energy_wca(energies(ENE%EXV))
 
+   if (flg_ele) call energy_ele_DH(energies(ENE%ELE))
+
    energies(ENE%TOTAL) = sum(energies(1:ENE%MAX))
 
 end subroutine energy
diff --git a/src/energy_ele_DH.F90 b/src/energy_ele_DH.F90
@@ -0,0 +1,36 @@
+subroutine energy_ele_DH(Eele)
+
+   use const
+   use pbc, only : pbc_vec_d
+   use var_state, only : xyz, lambdaD
+   use var_potential, only : ele_mp, nele, ele_cutoff, ele_coef
+
+   implicit none
+
+   real(PREC), intent(out) :: Eele
+
+   integer :: imp1, imp2, iele
+   real(PREC) :: dist
+   real(PREC) :: e_ele, rcdist
+
+   rcdist = 1.0_PREC / lambdaD
+   e_ele = 0.0_PREC
+
+   !$omp parallel do private(imp1,imp2,dist) reduction(+:e_ele)
+   do iele = 1, nele
+
+      imp1 = ele_mp(1, iele)
+      imp2 = ele_mp(2, iele)
+
+      dist = norm2(pbc_vec_d(xyz(:, ele_mp(1, iele)), xyz(:, ele_mp(2, iele))))
+
+      if (dist > ele_cutoff) cycle
+
+      e_ele = e_ele + ele_coef /dist*exp(-dist*rcdist)
+
+   end do
+   !$omp end parallel do
+
+   Eele = e_ele
+
+end subroutine energy_ele_DH
diff --git a/src/force.F90 b/src/force.F90
@@ -4,6 +4,7 @@ subroutine force()
    use const
    use const_idx, only : ENE
    use var_state, only : nthreads, forces
+   use var_potential, only : flg_ele
    use var_top, only : nmp
 
    implicit none
@@ -27,6 +28,7 @@ subroutine force()
    call force_angl(forces_t(1,1,tn))
    call force_bp(forces_t(1,1,tn))
    call force_wca(forces_t(1,1,tn))
+   if (flg_ele) call force_ele_DH(forces_t(1,1,tn))
 !$omp end parallel
 
    forces(:,:) = forces_t(:,:,0)

diff --git a/src/force_ele_DH.F90 b/src/force_ele_DH.F90
@@ -0,0 +1,55 @@
+subroutine force_ele_DH(forces)
+
+  use const, only : PREC
+  use pbc, only : pbc_vec_d
+  use var_top, only : nmp
+  use var_state, only : xyz, lambdaD
+  use var_potential, only : ele_coef, ele_cutoff, nele, ele_mp
+
+  implicit none
+
+  real(PREC), intent(inout) :: forces(3, nmp)
+
+  integer :: iele, imp1, imp2
+  real(PREC) :: dist1, dist2, rdist1
+  real(PREC) :: dvdw_dr, rcdist, cutoff2
+  real(PREC) :: v21(3), for(3)
+  !character(ARRAY_MSG_ERROR) :: error_message
+
+  ! --------------------------------------------------------------------
+
+  cutoff2 = ele_cutoff ** 2
+  rcdist = 1.0_PREC / lambdaD
+
+  !$omp do private(imp1,imp2,v21,dist2,dist1,rdist1,dvdw_dr,for)
+  do iele = 1, nele
+     imp1 = ele_mp(1, iele)
+     imp2 = ele_mp(2, iele)
+
+     v21(:) = pbc_vec_d(xyz(:,imp2), xyz(:,imp1))
+
+     dist2 = dot_product(v21,v21)
+
+     if(dist2 > cutoff2) cycle
+
+     ! -----------------------------------------------------------------
+     dist1 = sqrt(dist2)
+     rdist1 = 1.0 / dist1
+
+     dvdw_dr = ele_coef * rdist1 * rdist1 &
+              * (rdist1 + rcdist) * exp(-dist1 * rcdist)
+
+
+     !if(dvdw_dr > DE_MAX) then
+     !   write(error_message,*) 'force_ele_DH > DE_MAX', istep, grep, imp1, imp2, dist1, dvdw_dr, DE_MAX
+     !   write(6, '(a)') trim(error_message)
+     !   dvdw_dr = DE_MAX
+     !end if
+
+     for(1:3) = dvdw_dr * v21(1:3)
+     forces(1:3, imp1) = forces(1:3, imp1) - for(1:3)
+     forces(1:3, imp2) = forces(1:3, imp2) + for(1:3)
+  end do
+  !$omp end do nowait
+
+end subroutine force_ele_DH
diff --git a/src/job_dcd.F90 b/src/job_dcd.F90
@@ -38,7 +38,7 @@ subroutine job_dcd()
    allocate(xyz(3, nmp))
 
    open(hdl_out, file = cfile_out, status = 'replace', action = 'write', form='formatted')
-   write(hdl_out, '(a)') '#(1)nframe  (2)Etotal  (3)Ebond   (4)Eangl   (5)Ebp   (6)Eexv'
+   write(hdl_out, '(a)') '#(1)nframe (2)T   (3)Ekin       (4)Epot       (5)Ebond      (6)Eangl      (7)Ebp        (8)Eexv       (9)Eele'
 
    nframe = 0
    do
@@ -48,7 +48,7 @@ subroutine job_dcd()
 
       call energy()
 
-      write(hdl_out, '(i10, 5(1x,g15.8))') nframe, (energies(i), i=0,ENE%MAX)
+      write(hdl_out, '(i10, 1x, f6.2, 7(1x,g13.6))') nframe, 0.0, 0.0, (energies(i), i=0,ENE%MAX)
 
       !! Debug
       !if (nframe == 10) then

diff --git a/src/job_md.F90 b/src/job_md.F90
@@ -14,7 +14,7 @@ subroutine job_md()
                          flg_variable_box, variable_box_step, variable_box_change, &
                          opt_anneal, nanneal, anneal_tempK, anneal_step, &
                          istep, ianneal, istep_anneal_next
-   use var_potential, only : wca_nl_cut2, wca_sigma, bp_nl_cut2, bp_cutoff
+   use var_potential, only : wca_nl_cut2, wca_sigma, bp_nl_cut2, bp_cutoff, ele_cutoff, ele_nl_cut2
    use var_io, only : flg_progress, step_progress, hdl_dcd, hdl_out, cfile_prefix, cfile_out, cfile_pdb_ini, cfile_xyz_ini
    use dcd, only : file_dcd, DCD_OPEN_MODE
 
@@ -136,6 +136,7 @@ subroutine job_md()
    ! Neighbor list
    wca_nl_cut2 = (wca_sigma + nl_margin) ** 2
    bp_nl_cut2 = (bp_cutoff + nl_margin) ** 2
+   ele_nl_cut2 = (ele_cutoff + nl_margin) ** 2
    call neighbor_list()
    xyz_move(:,:) = 0.0e0_PREC
 
@@ -178,20 +179,20 @@ subroutine job_md()
 
    ! Open .out file
    open(hdl_out, file = cfile_out, status = 'replace', action = 'write', form='formatted')
-   write(hdl_out, '(a)') '#(1)nframe (2)T   (3)Ekin       (4)Epot       (5)Ebond      (6)Eangl      (7)Ebp        (8)Eexv'
-                         !1234567890 123456 1234567890123 1234567890123 1234567890123 1234567890123 1234567890123 1234567890123
+   write(hdl_out, '(a)') '#(1)nframe (2)T   (3)Ekin       (4)Epot       (5)Ebond      (6)Eangl      (7)Ebp        (8)Eexv       (9)Eele'
+                         !1234567890 123456 1234567890123 1234567890123 1234567890123 1234567890123 1234567890123 1234567890123 1234567890123
 
    ! Output initial structure
    if (restarted) then
       ! Only STDOUT if restarted. No DCD output.
       write(6, '(a)') '##### Energies at the beginning'
-      write(6, '(a)') '#(1)nframe (2)T   (3)Ekin       (4)Epot       (5)Ebond      (6)Eangl      (7)Ebp        (8)Eexv'
-      write(6, '(i10, 1x, f6.2, 6(1x,g13.6))') istep, tempK, Ekinetic, (energies(i), i=0,ENE%MAX)
+      write(6, '(a)') '#(1)nframe (2)T   (3)Ekin       (4)Epot       (5)Ebond      (6)Eangl      (7)Ebp        (8)Eexv       (9)Eele'
+      write(6, '(i10, 1x, f6.2, 7(1x,g13.6))') istep, tempK, Ekinetic, (energies(i), i=0,ENE%MAX)
       write(6, *)
 
    else
       ! At istep = 0 (not restarted), write both .out and DCD
-      write(hdl_out, '(i10, 1x, f6.2, 6(1x,g13.6))') istep, tempK, Ekinetic, (energies(i), i=0,ENE%MAX)
+      write(hdl_out, '(i10, 1x, f6.2, 7(1x,g13.6))') istep, tempK, Ekinetic, (energies(i), i=0,ENE%MAX)
       call fdcd%write_onestep(nmp, xyz)
    endif
 
@@ -263,7 +264,7 @@ subroutine job_md()
       if (mod(istep, nstep_save) == 0) then
          call energy()
          call energy_kinetic()
-         write(hdl_out, '(i10, 1x, f6.2, 6(1x,g13.6))') istep, tempK, Ekinetic, (energies(i), i=0,ENE%MAX)
+         write(hdl_out, '(i10, 1x, f6.2, 7(1x,g13.6))') istep, tempK, Ekinetic, (energies(i), i=0,ENE%MAX)
          call fdcd%write_onestep(nmp, xyz)
       endif
 

diff --git a/src/list_ele.F90 b/src/list_ele.F90
@@ -0,0 +1,64 @@
+subroutine list_ele()
+
+   use const
+   use var_top, only : nmp_chain, imp_chain, nchains, has_charge
+   use var_potential, only : nele, ele_mp
+
+   implicit none
+
+   integer :: ichain, jchain
+   integer :: n, i, j, j_start, imp, jmp
+   integer :: iele
+
+   if (allocated(ele_mp)) then
+      deallocate(ele_mp)
+   endif
+
+   do n = 1, 2
+
+      iele = 0
+
+      do ichain = 1, nchains 
+
+         do jchain = ichain, nchains 
+
+            do i = 1, nmp_chain(ichain)
+
+               imp = imp_chain(i, ichain)
+
+               if (.not. has_charge(imp)) cycle
+
+               if (ichain == jchain) then
+                  j_start = i + 3
+               else
+                  j_start = 1
+               endif
+
+               do j = j_start, nmp_chain(jchain)
+
+                  jmp = imp_chain(j, jchain)
+
+                  if (.not. has_charge(jmp)) cycle
+
+                  iele = iele + 1
+
+                  if (n == 2) then
+                     ele_mp(1,iele) = imp
+                     ele_mp(2,iele) = jmp
+                  endif
+
+               enddo
+            enddo
+
+         enddo
+      enddo
+
+      if (n == 1) then
+         nele = iele
+         allocate(ele_mp(2, nele))
+      endif
+   enddo
+
+   write(*,*) '#nele: ', nele
+
+end subroutine list_ele