Change FF file format to TOML

rna_cg1.ff

@@ -1,15 +1,28 @@
-bond_k  15.0
-bond_r0  5.9
-angl_k  10.0
-angl_t0  2.618
-Ubp_cutoff  18.0
-Ubp0  -1.66666666667
-Ubp_bond_k  3.0
-Ubp_bond_r  13.8
-Ubp_angl_k  1.5
-Ubp_angl_theta1  1.8326
-Ubp_angl_theta2  0.9425
-Ubp_dihd_k  0.5
-Ubp_dihd_phi1  1.8326
-Ubp_dihd_phi2  1.1345
-Ubp_min_loop 4
+[potential]
+    [potential.bond]
+        k = 15.0
+        r0 = 5.9
+
+    [potential.angle]
+        k = 10.0
+        a0 = 2.618
+
+    [potential.basepair]
+        min_loop = 4
+        cutoff = 18.0
+        U0 = -1.66666666667
+
+        bond_k = 3.0
+        bond_r = 13.8
+
+        angl_k = 1.5
+        angl_theta1 = 1.8326
+        angl_theta2 = 0.9425
+
+        dihd_k = 0.5
+        dihd_phi1 = 1.8326
+        dihd_phi2 = 1.1345
+
+    [potential.wca]
+        sigma = 10.0
+        epsilon = 2.0

rna_cg1_para2.ff

@@ -1,15 +1,30 @@
-bond_k  15.0
-bond_r0  5.9
-angl_k  10.0
-angl_t0  2.618
-Ubp_cutoff  18.0
-Ubp0  -2.45
-Ubp_bond_k  3.0
-Ubp_bond_r  13.8
-Ubp_angl_k  3.2
-Ubp_angl_theta1  1.8326
-Ubp_angl_theta2  0.9425
-Ubp_dihd_k  1.3
-Ubp_dihd_phi1  1.8326
-Ubp_dihd_phi2  1.1345
-Ubp_min_loop 3
+title = "sis force field refined parameter 2"
+
+[potential]
+    [potential.bond]
+        k = 15.0
+        r0 = 5.9
+
+    [potential.angle]
+        k = 10.0
+        a0 = 2.618
+
+    [potential.basepair]
+        min_loop = 3
+        cutoff = 18.0
+        U0 = -2.45
+
+        bond_k = 3.0
+        bond_r = 13.8
+
+        angl_k = 3.2
+        angl_theta1 = 1.8326
+        angl_theta2 = 0.9425
+
+        dihd_k = 1.3
+        dihd_phi1 = 1.8326
+        dihd_phi2 = 1.1345
+
+    [potential.wca]
+        sigma = 10.0
+        epsilon = 2.0

src/energy_bp.F90

@@ -26,29 +26,29 @@ subroutine energy_bp(Ebp)
 
       d = mp_distance(imp, jmp)
 
-      if (d >= Ubp_cutoff) cycle
+      if (d >= bp_cutoff) cycle
 
-      u = Ubp_bond_k * (d - Ubp_bond_r)**2
+      u = bp_bond_k * (d - bp_bond_r)**2
 
       theta = mp_angle(imp, jmp, jmp-1)
-      u = u + Ubp_angl_k * (theta - Ubp_angl_theta1)**2
+      u = u + bp_angl_k * (theta - bp_angl_theta1)**2
 
       theta = mp_angle(imp-1, imp, jmp)
-      u = u + Ubp_angl_k * (theta - Ubp_angl_theta1)**2
+      u = u + bp_angl_k * (theta - bp_angl_theta1)**2
 
       theta = mp_angle(imp, jmp, jmp+1)
-      u = u + Ubp_angl_k * (theta - Ubp_angl_theta2)**2
+      u = u + bp_angl_k * (theta - bp_angl_theta2)**2
 
       theta = mp_angle(imp+1, imp, jmp)
-      u = u + Ubp_angl_k * (theta - Ubp_angl_theta2)**2
+      u = u + bp_angl_k * (theta - bp_angl_theta2)**2
 
       phi = mp_dihedral(imp-1, imp, jmp, jmp-1)
-      u = u + Ubp_dihd_k * (1.0 + cos(phi + Ubp_dihd_phi1))
+      u = u + bp_dihd_k * (1.0 + cos(phi + bp_dihd_phi1))
 
       phi = mp_dihedral(imp+1, imp, jmp, jmp+1)
-      u = u + Ubp_dihd_k * (1.0 + cos(phi + Ubp_dihd_phi2))
+      u = u + bp_dihd_k * (1.0 + cos(phi + bp_dihd_phi2))
 
-      e_bp(ibp) = nhb * Ubp0 * exp(-u)
+      e_bp(ibp) = nhb * bp_U0 * exp(-u)
 
    enddo
    !$omp end parallel do

src/force_bp.F90

@@ -42,18 +42,18 @@ subroutine force_bp()
       d1212 = dot_product(v12,v12)
       a12 = sqrt(d1212)
 
-      if (a12 >= Ubp_cutoff) cycle
+      if (a12 >= bp_cutoff) cycle
 
       imp3 = imp1 - 1
       imp4 = imp2 - 1
       imp5 = imp1 + 1
       imp6 = imp2 + 1
 
       !===== Distance =====
-      d = a12 - Ubp_bond_r
+      d = a12 - bp_bond_r
 
-      u = Ubp_bond_k * d**2
-      f_i(:) = (2.0e0_PREC * Ubp_bond_k * d / a12) * v12(:)
+      u = bp_bond_k * d**2
+      f_i(:) = (2.0e0_PREC * bp_bond_k * d / a12) * v12(:)
       f_bp(:, 1, ibp) = + f_i(:)
       f_bp(:, 2, ibp) = - f_i(:)
 
@@ -77,9 +77,9 @@ subroutine force_bp()
 
       !===== Angle of 3-1=2 (imp-1 -- imp -- jmp) =====
       cosine = d1213 / (sqrt(d1313) * a12)
-      d = acos(cosine) - Ubp_angl_theta1
-      pre = 2.0e0_PREC * Ubp_angl_k * d / sqrt(d1313*d1212 - d1213**2)
-      u = u + Ubp_angl_k * d**2
+      d = acos(cosine) - bp_angl_theta1
+      pre = 2.0e0_PREC * bp_angl_k * d / sqrt(d1313*d1212 - d1213**2)
+      u = u + bp_angl_k * d**2
 
       f_i(:) = pre * (v12(:) - (d1213 / d1313 * v13(:)))
       f_k(:) = pre * (v13(:) - (d1213over1212 * v12(:)))
@@ -90,9 +90,9 @@ subroutine force_bp()
 
       !===== Angle of 1=2-4 (imp -- jmp -- jmp-1) =====
       cosine = d1242 / (a12 * sqrt(d4242))
-      d = acos(cosine) - Ubp_angl_theta1
-      u = u + Ubp_angl_k * d**2
-      pre = 2.0e0_PREC * Ubp_angl_k * d / sqrt(d1212*d4242 - d1242**2)
+      d = acos(cosine) - bp_angl_theta1
+      u = u + bp_angl_k * d**2
+      pre = 2.0e0_PREC * bp_angl_k * d / sqrt(d1212*d4242 - d1242**2)
 
       f_i(:) = - pre * (v42(:) - (d1242over1212 * v12(:)))
       f_k(:) = - pre * (v12(:) - (d1242 / d4242 * v42(:)))
@@ -103,9 +103,9 @@ subroutine force_bp()
 
       !===== Angle of 5-1=2 (imp+1 -- imp -- jmp) =====
       cosine = d1215 / (sqrt(d1515) * a12)
-      d = acos(cosine) - Ubp_angl_theta2
-      u = u + Ubp_angl_k * d**2
-      pre = 2.0e0_PREC * Ubp_angl_k * d / sqrt(d1515*d1212 - d1215**2)
+      d = acos(cosine) - bp_angl_theta2
+      u = u + bp_angl_k * d**2
+      pre = 2.0e0_PREC * bp_angl_k * d / sqrt(d1515*d1212 - d1215**2)
 
       f_i(:) = pre * (v12(:) - (d1215 / d1515 * v15(:)))
       f_k(:) = pre * (v15(:) - (d1215over1212 * v12(:)))
@@ -116,9 +116,9 @@ subroutine force_bp()
 
       !===== Angle of 1=2-6 (imp -- jmp -- jmp+1) =====
       cosine = d1262 / (a12 * sqrt(d6262))
-      d = acos(cosine) - Ubp_angl_theta2
-      u = u + Ubp_angl_k * d**2
-      pre = 2.0e0_PREC * Ubp_angl_k * d / sqrt(d1212*d6262 - d1262**2)
+      d = acos(cosine) - bp_angl_theta2
+      u = u + bp_angl_k * d**2
+      pre = 2.0e0_PREC * bp_angl_k * d / sqrt(d1212*d6262 - d1262**2)
 
       f_i(:) = - pre * (v62(:) - (d1262over1212 * v12(:)))
       f_k(:) = - pre * (v12(:) - (d1262 / d6262 * v62(:)))
@@ -136,10 +136,10 @@ subroutine force_bp()
       n(3) = v12(1)*v13(2) - v12(2)*v13(1)
 
       dih = atan2(dot_product(v42,n)*a12, dot_product(m,n))
-      d = dih + Ubp_dihd_phi1
-      u = u + Ubp_dihd_k * (1.0 + cos(d))
+      d = dih + bp_dihd_phi1
+      u = u + bp_dihd_k * (1.0 + cos(d))
 
-      pre = -Ubp_dihd_k * sin(d) * a12
+      pre = -bp_dihd_k * sin(d) * a12
       f_i(:) = + pre / dot_product(m, m) * m(:)
       f_l(:) = - pre / dot_product(n, n) * n(:)
 
@@ -159,10 +159,10 @@ subroutine force_bp()
       n(3) = v12(1)*v15(2) - v12(2)*v15(1)
 
       dih = atan2(dot_product(v62,n)*a12, dot_product(m,n))
-      d = dih + Ubp_dihd_phi2
-      u = u + Ubp_dihd_k * (1.0 + cos(d))
+      d = dih + bp_dihd_phi2
+      u = u + bp_dihd_k * (1.0 + cos(d))
 
-      pre = -Ubp_dihd_k * sin(d) * a12
+      pre = -bp_dihd_k * sin(d) * a12
       f_i(:) = + pre / dot_product(m, m) * m(:)
       f_l(:) = - pre / dot_product(n, n) * n(:)
 
@@ -174,7 +174,7 @@ subroutine force_bp()
       f_bp(:, 5, ibp) = f_bp(:, 5, ibp) + f_l(:)
 
       !===== Total =====
-      f_bp(:, :, ibp) = bp_mp(3, ibp) * Ubp0 * exp(-u) * f_bp(:, :, ibp)
+      f_bp(:, :, ibp) = bp_mp(3, ibp) * bp_U0 * exp(-u) * f_bp(:, :, ibp)
 
    enddo
    !$omp end parallel do

src/list_bp.F90

@@ -3,7 +3,7 @@ subroutine list_bp()
    use const
    use const_idx, only : SEQT
    use var_top, only : nmp_chain, seq, imp_chain, nchains
-   use var_potential, only : nbp, bp_mp, Ubp_min_loop
+   use var_potential, only : nbp, bp_mp, bp_min_loop
 
    implicit none
 
@@ -24,7 +24,7 @@ subroutine list_bp()
                imp = imp_chain(i, ichain)
 
                if (ichain == jchain) then
-                  j_start = i + Ubp_min_loop + 1
+                  j_start = i + bp_min_loop + 1
                else
                   j_start = 2
                endif