Skip to content

Commit

Permalink
Change BP parameter format to be fully nucleotide-pair dependent.
Browse files Browse the repository at this point in the history
  • Loading branch information
@naotohori
naotohori committed Jun 2, 2022
1 parent 2adda05 commit b663cc4
Show file tree
Hide file tree
Showing 14 changed files with 432 additions and 275 deletions.
43 changes: 36 additions & 7 deletions check/rna_cg1_para3.ff
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,20 +10,49 @@ title = "sis force field refined parameter 2"
a0 = 2.618

[potential.basepair]
min_loop = 3
cutoff = 18.0
U0_GC = -8.40
U0_AU = -4.20
U0_GU = -2.80

[potential.basepair.GC]
U0 = -8.40
bond_k = 3.0
bond_r = 13.8

angl_k = 3.2
angl_k1 = 3.2
angl_k2 = 3.2
angl_theta1 = 1.8326
angl_theta2 = 0.9425

dihd_k = 1.3
dihd_k1 = 1.3
dihd_k2 = 1.3
dihd_phi1 = 1.8326
dihd_phi2 = 1.1345

[potential.basepair.AU]
U0 = -4.20
bond_k = 3.0
bond_r = 13.8

angl_k1 = 3.2
angl_k2 = 3.2
angl_theta1 = 1.8326
angl_theta2 = 0.9425

dihd_k1 = 1.3
dihd_k2 = 1.3
dihd_phi1 = 1.8326
dihd_phi2 = 1.1345

[potential.basepair.GU]
U0 = -2.80
bond_k = 3.0
bond_r = 13.8

angl_k1 = 3.2
angl_k2 = 3.2
angl_theta1 = 1.8326
angl_theta2 = 0.9425

dihd_k1 = 1.3
dihd_k2 = 1.3
dihd_phi1 = 1.8326
dihd_phi2 = 1.1345

Expand Down
47 changes: 36 additions & 11 deletions para/rna_NHT22.ff
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,20 +11,45 @@ title = 'NHG22'

[potential.basepair]
cutoff_energy = -0.01
U0_GC = -5.00000
U0_AU = -3.33333
U0_GU = -3.33333

bond_k = 3.0
bond_r = 13.8
[potential.basepair.GC]
U0 = -5.00000
bond_k = 3.0
bond_r = 13.8
angl_k1 = 1.5
angl_k2 = 1.5
angl_theta1 = 1.8326
angl_theta2 = 0.9425
dihd_k1 = 0.5
dihd_k2 = 0.5
dihd_phi1 = 1.8326
dihd_phi2 = 1.1345

angl_k = 1.5
angl_theta1 = 1.8326
angl_theta2 = 0.9425
[potential.basepair.AU]
U0 = -3.33333
bond_k = 3.0
bond_r = 13.8
angl_k1 = 1.5
angl_k2 = 1.5
angl_theta1 = 1.8326
angl_theta2 = 0.9425
dihd_k1 = 0.5
dihd_k2 = 0.5
dihd_phi1 = 1.8326
dihd_phi2 = 1.1345

dihd_k = 0.5
dihd_phi1 = 1.8326
dihd_phi2 = 1.1345
[potential.basepair.GU]
U0 = -3.33333
bond_k = 3.0
bond_r = 13.8
angl_k1 = 1.5
angl_k2 = 1.5
angl_theta1 = 1.8326
angl_theta2 = 0.9425
dihd_k1 = 0.5
dihd_k2 = 0.5
dihd_phi1 = 1.8326
dihd_phi2 = 1.1345

[potential.wca]
sigma = 10.0
Expand Down
10 changes: 6 additions & 4 deletions src/const_idx.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,14 +43,16 @@ module const_idx
type(integrator_types), parameter :: INTGRT = integrator_types(0,1,1)

type bp_types
integer :: UNDEF ! 0
integer :: GC_WCF ! 1
integer :: AU_WCF ! 2
integer :: GU_WBL ! 3
integer :: UNDEF ! 0
integer :: GC ! 1
integer :: AU ! 2
integer :: GU ! 3
integer :: MAX
endtype bp_types
type(bp_types), parameter :: BPT = bp_types(0,1,2,3,3)

character(len=*), dimension(3), parameter :: BPTYPE_CHAR = (/'GC', 'AU', 'GU'/)

type rst_block
integer :: STEP
integer :: ANNEAL
Expand Down
28 changes: 15 additions & 13 deletions src/energy_bp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,51 +4,53 @@ subroutine energy_bp(Ebp)
use const_phys, only : ZERO_JUDGE
use pbc, only : pbc_vec_d
use var_state, only : xyz, kT
use var_potential
use var_potential, only : bp_paras, nbp, bp_mp, bp_type2nhb, basepair_parameters
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

real(PREC), intent(inout) :: Ebp

integer :: ibp, imp, jmp, nhb
type(basepair_parameters) :: bpp
real(PREC) :: u
real(PREC) :: d, theta, phi
real(PREC) :: e_bp(nbp)

e_bp(:) = 0.0e0_PREC

!$omp parallel do private(imp,jmp,d,u,theta,phi)
!$omp parallel do private(imp, jmp, d, u, theta, phi, bpp)
do ibp = 1, nbp

imp = bp_mp(1, ibp)
jmp = bp_mp(2, ibp)
bpp = bp_paras(bp_mp(3, ibp))

d = norm2(pbc_vec_d(xyz(:,imp), xyz(:, jmp))) - bp_bond_r
d = norm2(pbc_vec_d(xyz(:,imp), xyz(:, jmp))) - bpp%bond_r

if (abs(d) > bp_cutoff_ddist) cycle
u = bp_bond_k * d**2
if (abs(d) > bpp%cutoff_ddist) cycle

u = bpp%bond_k * d**2

theta = mp_angle(imp, jmp, jmp-1)
u = u + bp_angl_k * (theta - bp_angl_theta1)**2
u = u + bpp%angl_k1 * (theta - bpp%angl_theta1)**2

theta = mp_angle(imp-1, imp, jmp)
u = u + bp_angl_k * (theta - bp_angl_theta1)**2
u = u + bpp%angl_k1 * (theta - bpp%angl_theta1)**2

theta = mp_angle(imp, jmp, jmp+1)
u = u + bp_angl_k * (theta - bp_angl_theta2)**2
u = u + bpp%angl_k2 * (theta - bpp%angl_theta2)**2

theta = mp_angle(imp+1, imp, jmp)
u = u + bp_angl_k * (theta - bp_angl_theta2)**2
u = u + bpp%angl_k2 * (theta - bpp%angl_theta2)**2

phi = mp_dihedral(imp-1, imp, jmp, jmp-1)
u = u + bp_dihd_k * (1.0 + cos(phi + bp_dihd_phi1))
u = u + bpp%dihd_k1 * (1.0_PREC + cos(phi + bpp%dihd_phi1))

phi = mp_dihedral(imp+1, imp, jmp, jmp+1)
u = u + bp_dihd_k * (1.0 + cos(phi + bp_dihd_phi2))
u = u + bpp%dihd_k2 * (1.0_PREC + cos(phi + bpp%dihd_phi2))

e_bp(ibp) = bp_U0(ibp) * exp(-u)
e_bp(ibp) = bpp%U0 * exp(-u)

enddo
!$omp end parallel do
Expand Down
28 changes: 14 additions & 14 deletions src/energy_bp_limit.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,7 @@ subroutine energy_bp_limit(Ebp)
use pbc, only : pbc_vec_d
use var_top, only : nmp
use var_state, only : xyz, kT, bp_status, ene_bp, flg_bp_energy
use var_potential, only : max_bp_per_nt, bp_cutoff_ddist, bp_cutoff_ene, bp_bond_k, bp_bond_r, &
bp_angl_k, bp_angl_theta1, bp_angl_theta2, bp_dihd_k, bp_dihd_phi1, bp_dihd_phi2, &
nbp, bp_mp, bp_U0
use var_potential, only : max_bp_per_nt, bp_cutoff_ene, nbp, bp_mp, bp_paras, basepair_parameters
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
Expand All @@ -18,6 +16,7 @@ subroutine energy_bp_limit(Ebp)
integer :: nt_delete, ibp_delete
integer :: imp, jmp
integer :: i_save, i_swap
type(basepair_parameters) :: bpp
real(PREC) :: u, beta, ratio
real(PREC) :: d, theta, phi
real(PREC) :: ene
Expand All @@ -38,37 +37,38 @@ subroutine energy_bp_limit(Ebp)

!$omp barrier

!$omp parallel do private(imp, jmp, d, u, theta, phi, ene)
!$omp parallel do private(imp, jmp, d, u, theta, phi, ene, bpp)
do ibp = 1, nbp

imp = bp_mp(1, ibp)
jmp = bp_mp(2, ibp)
bpp = bp_paras(bp_mp(3, ibp))

d = norm2(pbc_vec_d(xyz(:,imp), xyz(:, jmp))) - bp_bond_r
d = norm2(pbc_vec_d(xyz(:,imp), xyz(:, jmp))) - bpp%bond_r

if (abs(d) > bp_cutoff_ddist) cycle
if (abs(d) > bpp%cutoff_ddist) cycle

u = bp_bond_k * d**2
u = bpp%bond_k * d**2

theta = mp_angle(imp, jmp, jmp-1)
u = u + bp_angl_k * (theta - bp_angl_theta1)**2
u = u + bpp%angl_k1 * (theta - bpp%angl_theta1)**2

theta = mp_angle(imp-1, imp, jmp)
u = u + bp_angl_k * (theta - bp_angl_theta1)**2
u = u + bpp%angl_k1 * (theta - bpp%angl_theta1)**2

theta = mp_angle(imp, jmp, jmp+1)
u = u + bp_angl_k * (theta - bp_angl_theta2)**2
u = u + bpp%angl_k2 * (theta - bpp%angl_theta2)**2

theta = mp_angle(imp+1, imp, jmp)
u = u + bp_angl_k * (theta - bp_angl_theta2)**2
u = u + bpp%angl_k2 * (theta - bpp%angl_theta2)**2

phi = mp_dihedral(imp-1, imp, jmp, jmp-1)
u = u + bp_dihd_k * (1.0 + cos(phi + bp_dihd_phi1))
u = u + bpp%dihd_k1 * (1.0_PREC + cos(phi + bpp%dihd_phi1))

phi = mp_dihedral(imp+1, imp, jmp, jmp+1)
u = u + bp_dihd_k * (1.0 + cos(phi + bp_dihd_phi2))
u = u + bpp%dihd_k2 * (1.0_PREC + cos(phi + bpp%dihd_phi2))

ene = bp_U0(ibp) * exp(-u)
ene = bpp%U0 * exp(-u)

if (ene <= bp_cutoff_ene) then
ene_bp(ibp) = ene
Expand Down
Loading

0 comments on commit b663cc4

Please sign in to comment.