Skip to content

Commit

Permalink
ANNaMo: add the possibility of using different interaction parameters…
Browse files Browse the repository at this point in the history 
… for crossover beads
  • Loading branch information
@lorenzo-rovigatti
lorenzo-rovigatti committed Oct 18, 2024
1 parent 8ba1129 commit 8ad9ce2
Show file tree
Hide file tree
Showing 4 changed files with 82 additions and 72 deletions.
103 changes: 60 additions & 43 deletions contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,8 @@ void CGNucleicAcidsInteraction::get_settings(input_file &inp) {
getInputNumber(&inp, "DPS_rfene", &_rfene, 0);
getInputNumber(&inp, "DPS_Kfene", &_Kfene, 0);
getInputNumber(&inp, "DPS_WCA_sigma", &_WCA_sigma, 0);
getInputNumber(&inp, "DPS_WCA_sigma_crossover", &_WCA_sigma_crossover, 0);
getInputString(&inp, "DPS_crossover_file", _crossover_file, 0);

if(getInputNumber(&inp, "max_backbone_force", &_mbf_fmax, 0) == KEY_FOUND) {
_use_mbf = true;
Expand All @@ -82,6 +84,7 @@ void CGNucleicAcidsInteraction::init() {
_PS_sqr_rep_rcut = pow(2., 2. / _PS_n) * SQR(_WCA_sigma);
_WCA_sigma_unbonded = _WCA_sigma * (6.0 / _bead_size - _3b_sigma) / 2.0; // not disabled
_PS_sqr_rep_rcut_unbonded = pow(2., 2. / _PS_n) * SQR(_WCA_sigma_unbonded);
_PS_sqr_rep_rcut_crossover = pow(2., 2. / _PS_n) * SQR(_WCA_sigma_crossover);

OX_LOG(Logger::LOG_INFO, "CGNA: WCA sigma = %lf, WCA sigma unbonded = %lf", _WCA_sigma, _WCA_sigma_unbonded);

Expand Down Expand Up @@ -215,9 +218,20 @@ number CGNucleicAcidsInteraction::_fene(BaseParticle *p, BaseParticle *q, bool u

number CGNucleicAcidsInteraction::_WCA(BaseParticle *p, BaseParticle *q, bool update_forces) {
number sqr_r = _computed_r.norm();

number sigma = p->is_bonded(q) ? _WCA_sigma : _WCA_sigma_unbonded;
number sqr_rcut = p->is_bonded(q) ? _PS_sqr_rep_rcut : _PS_sqr_rep_rcut_unbonded;

number sigma = _WCA_sigma;
number sqr_rcut = _PS_sqr_rep_rcut;

if(!p->is_bonded(q)) {
if(crossovers[p->index] || crossovers[q->index]) {
sigma = _WCA_sigma_crossover;
sqr_rcut = _PS_sqr_rep_rcut_crossover;
}
else {
sigma = _WCA_sigma_unbonded;
sqr_rcut = _PS_sqr_rep_rcut_unbonded;
}
}

if(sqr_r > sqr_rcut) {
return (number) 0.;
Expand All @@ -227,8 +241,6 @@ number CGNucleicAcidsInteraction::_WCA(BaseParticle *p, BaseParticle *q, bool up
// this number is the module of the force over r, so we don't have to divide the distance vector for its module
number force_mod = 0;



// cut-off for all the repulsive interactions
if(sqr_r < sqr_rcut) {
number part = 1.;
Expand All @@ -241,13 +253,6 @@ number CGNucleicAcidsInteraction::_WCA(BaseParticle *p, BaseParticle *q, bool up
force_mod += 4. * _PS_n * part * (2. * part - 1.) / sqr_r;
}
}
// attraction
/*else if(int_type == 0 && _PS_alpha != 0.) {
energy += 0.5 * _PS_alpha * (cos(_PS_gamma * sqr_r + _PS_beta) - 1.0);
if(update_forces) {
force_mod += _PS_alpha * _PS_gamma * sin(_PS_gamma * sqr_r + _PS_beta);
}
}*/

if(update_forces) {
auto force = -_computed_r * force_mod;
Expand Down Expand Up @@ -371,7 +376,6 @@ number CGNucleicAcidsInteraction::_patchy_three_body(BaseParticle *p, PSBond &ne
energy += prefactor * curr_energy * other_energy;

if(update_forces) {
// if(new_bond.r_mod > _3b_sigma)
{
BaseParticle *other = new_bond.other;

Expand All @@ -393,7 +397,6 @@ number CGNucleicAcidsInteraction::_patchy_three_body(BaseParticle *p, PSBond &ne
_update_stress_tensor(p->pos + new_bond.r, -tmp_force);
}

// if(other_bond.r_mod > _3b_sigma)
{
BaseParticle *other = other_bond.other;

Expand Down Expand Up @@ -584,20 +587,7 @@ number CGNucleicAcidsInteraction::pair_nonbonded_WCA(BaseParticle *p, BasePartic
number energy = _WCA(p, q, update_forces);

return energy;
} // switch(N_int_centers()) {
// case 0:
// break;
// case 1: {
// _base_patches[0] = LR_vector(1, 0, 0);
// break;
// }
// default:
// throw oxDNAException("Unsupported number of patches %d\n", N_int_centers());
// }
// for(uint i = 0; i < N_int_centers(); i++) {
// _base_patches[i].normalize();
// _base_patches[i] *= deltaPM;
// }
}

number CGNucleicAcidsInteraction::pair_nonbonded_sticky(BaseParticle *p, BaseParticle *q, bool compute_r, bool update_forces) {
if(p->is_bonded(q)) {
Expand Down Expand Up @@ -638,7 +628,7 @@ void CGNucleicAcidsInteraction::_parse_interaction_matrix() {
const number _kB_ = 1.9872036;
number _mu = _t37_ / (_tC+273.15);
if(inter_matrix_file.state == ERROR) {
throw oxDNAException("Caught an error while opening the interaction matrix file '%s'", _interaction_matrix_file.c_str());
throw oxDNAException("Caught an error while opening the interaction matrix file '%s'", _interaction_matrix_file.c_str());
}

_interaction_matrix_size = _N_attractive_types + 1;
Expand All @@ -647,20 +637,20 @@ void CGNucleicAcidsInteraction::_parse_interaction_matrix() {
ofstream myfile;
myfile.open("beta_eps_matrix.dat");
for(int i = 1; i <= _N_attractive_types; i++) {
for(int j = 1; j <= _N_attractive_types; j++) {
number valueH;
number valueS;
std::string keyH = Utils::sformat("dH[%d][%d]", i, j);
std::string keyS = Utils::sformat("dS[%d][%d]", i, j);
if(getInputNumber(&inter_matrix_file, keyH.c_str(), &valueH, 0) == KEY_FOUND && getInputNumber(&inter_matrix_file, keyS.c_str(), &valueS, 0) == KEY_FOUND) {
number beta_dG = (_mu * valueH * 1000 / _t37_ - valueS)/_kB_;
number beta_eps = -(beta_dG + dS_mod) / alpha_mod;
if(beta_dG<bdG_threshold && abs(i - j) > 1) {
_3b_epsilon[i + _interaction_matrix_size * j] = _3b_epsilon[j + _interaction_matrix_size * i] = beta_eps;
myfile << "beta_eps[" << i << "][" << j << "]=" << beta_eps << "\n";
}
}
for(int j = 1; j <= _N_attractive_types; j++) {
number valueH;
number valueS;
std::string keyH = Utils::sformat("dH[%d][%d]", i, j);
std::string keyS = Utils::sformat("dS[%d][%d]", i, j);
if(getInputNumber(&inter_matrix_file, keyH.c_str(), &valueH, 0) == KEY_FOUND && getInputNumber(&inter_matrix_file, keyS.c_str(), &valueS, 0) == KEY_FOUND) {
number beta_dG = (_mu * valueH * 1000 / _t37_ - valueS)/_kB_;
number beta_eps = -(beta_dG + dS_mod) / alpha_mod;
if(beta_dG<bdG_threshold && abs(i - j) > 1) {
_3b_epsilon[i + _interaction_matrix_size * j] = _3b_epsilon[j + _interaction_matrix_size * i] = beta_eps;
myfile << "beta_eps[" << i << "][" << j << "]=" << beta_eps << "\n";
}
}
}
}
myfile.close();
}
Expand All @@ -670,6 +660,33 @@ void CGNucleicAcidsInteraction::read_topology(int *N_strands, std::vector<BasePa
unsigned int N_from_conf = particles.size();
BaseInteraction::read_topology(N_strands, particles);

crossovers = std::vector<int>(N_from_conf, 0);
if(_crossover_file.size() > 0) {
std::ifstream c_file(_crossover_file.c_str());
if(!c_file) {
throw oxDNAException("The crossover file '%s' does not exist or is unreadable", _crossover_file.c_str());
}

int crossover_read = 0;
while(c_file.good()) {
std::getline(c_file, line);
auto spl = Utils::split(line);
if(spl.size() == 1) {
unsigned int idx = std::atoi(spl[0].c_str());
if(idx >= N_from_conf) {
throw oxDNAException("The index %u found in the crossover file exceeds the number of particles %u", idx, N_from_conf);
}
crossovers[idx] = 1;
crossover_read++;
}
else if(spl.size() > 1) {
throw oxDNAException("Malformed crossover file: each line should have no more than one column");
}
}
OX_LOG(Logger::LOG_INFO, "Read %d crossovers from the file", crossover_read);
c_file.close();
}

std::ifstream topology;
topology.open(_topology_filename, std::ios::in);
if(!topology.good()) {
Expand Down Expand Up @@ -767,7 +784,7 @@ void CGNucleicAcidsInteraction::read_topology(int *N_strands, std::vector<BasePa
_N_chains = *N_strands;

_parse_interaction_matrix();

return;
}

Expand Down
6 changes: 4 additions & 2 deletions contrib/tostiguerra/src/CGNucleicAcidsInteraction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,8 @@ class CGNucleicAcidsInteraction: public BaseInteraction {
number _Kfene = 15.;
number _rfene = 1.5;
number _sqr_rfene;
number _WCA_sigma = 1.0, _WCA_sigma_unbonded;
number _PS_sqr_rep_rcut, _PS_sqr_rep_rcut_unbonded;
number _WCA_sigma = 1.0, _WCA_sigma_unbonded, _WCA_sigma_crossover = 0.0;
number _PS_sqr_rep_rcut, _PS_sqr_rep_rcut_unbonded, _PS_sqr_rep_rcut_crossover;
number _tC = 37.0;
number dS_mod = 1.0;
number alpha_mod = 1.0;
Expand All @@ -65,6 +65,7 @@ class CGNucleicAcidsInteraction: public BaseInteraction {
int _N_attractive_types = 0;
int _interaction_matrix_size = 0;
std::string _interaction_matrix_file;
std::string _crossover_file;

/// three-body potential stuff
number _3b_rcut = -1.;
Expand Down Expand Up @@ -98,6 +99,7 @@ class CGNucleicAcidsInteraction: public BaseInteraction {
number _mbf_A = 0.0;
number _mbf_B = 0.0;

std::vector<int> crossovers;
std::map<int, std::set<PSBond, PSBondCompare> > _bonds;

int _N_chains = -1;
Expand Down
44 changes: 17 additions & 27 deletions contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,10 +22,12 @@ __constant__ int MD_n[1];
__constant__ int MD_interaction_matrix_size[1];
__constant__ float MD_sqr_rep_rcut[1];
__constant__ float MD_sqr_rep_rcut_unbonded[1];
__constant__ float MD_sqr_rep_rcut_crossover[1];
__constant__ float MD_sqr_rfene[1];
__constant__ float MD_Kfene[1];
__constant__ float MD_WCA_sigma[1];
__constant__ float MD_WCA_sigma_unbonded[1];
__constant__ float MD_WCA_sigma_crossover[1];
__constant__ float MD_deltaPatchMon[1];
__constant__ float MD_sqr_rcut[1];
__constant__ float MD_alpha[1];
Expand Down Expand Up @@ -205,7 +207,6 @@ __device__ void _patchy_three_body(CUDA_FS_bond_list &bond_list, c_number4 &F, c
// the factor 2 takes into account the fact that the pair energy is always counted twice
F.w += 2.f * prefactor * curr_energy * other_energy;

// if(new_bond.r_mod > _3b_sigma)
{
c_number factor = -prefactor * other_energy;
c_number4 tmp_force = factor * b1.force;
Expand All @@ -220,7 +221,6 @@ __device__ void _patchy_three_body(CUDA_FS_bond_list &bond_list, c_number4 &F, c
// _update_stress_tensor(p->pos + bi.r, -tmp_force);
}

// if(other_bond.r_mod > _3b_sigma)
{
c_number factor = -prefactor * curr_energy;
c_number4 tmp_force = factor * b2.force;
Expand All @@ -234,26 +234,6 @@ __device__ void _patchy_three_body(CUDA_FS_bond_list &bond_list, c_number4 &F, c
// _update_stress_tensor(p->pos, tmp_force);
// _update_stress_tensor(p->pos + bj.r, -tmp_force);
}

// if(curr_energy != MD_3b_epsilon[0]) {
// c_number factor = MD_3b_prefactor[0] * other_energy;
// c_number4 force = b1.force * factor;
// force.w = 0.f;
//
// _update_stress_tensor<false>(p_st, b1.r, force);
// F += force;
// LR_atomicAddXYZ(forces + b1.q, -force);
// }
//
// if(other_energy != MD_3b_epsilon[0]) {
// c_number factor = MD_3b_prefactor[0] * curr_energy;
// c_number4 force = b2.force * factor;
// force.w = 0.f;
//
// _update_stress_tensor<false>(p_st, b2.r, force);
// F += force;
// LR_atomicAddXYZ(forces + b2.q, -force);
// }
}
}
}
Expand All @@ -273,7 +253,7 @@ __device__ void _flexibility_three_body(c_number4 &ppos, c_number4 &n1_pos, c_nu
c_number force_mod_n2 = i_pn1_pn2 + cost_n2;

c_number energy, force_factor;
if(MD_semiflexibility_3b_exp_sigma[0] > 0.0) {
if(MD_semiflexibility_3b_exp_sigma[0] > 0.0f) {
c_number arg = (1.f - cost) / MD_semiflexibility_3b_exp_sigma[0];
c_number exp_factor = expf(-SQR(arg));
energy = -MD_semiflexibility_3b_k[0] * (exp_factor - 1.f);
Expand Down Expand Up @@ -361,8 +341,9 @@ __device__ bool _sticky_interaction(int p_btype, int q_btype) {
return true;
}

__global__ void ps_forces(c_number4 *poss, GPU_quat *orientations, c_number4 *forces, c_number4 *torques, c_number4 *three_body_forces, c_number4 *three_body_torques, int *matrix_neighs, int *number_neighs, bool update_st,
cudaTextureObject_t tex_eps, CUDAStressTensor *st, CUDABox *box) {
__global__ void ps_forces(c_number4 *poss, GPU_quat *orientations, c_number4 *forces, c_number4 *torques, c_number4 *three_body_forces,
c_number4 *three_body_torques, int *matrix_neighs, int *number_neighs, int *crossovers, bool update_st,
cudaTextureObject_t tex_eps, CUDAStressTensor *st, CUDABox *box) {
if(IND >= MD_N[0]) return;

c_number4 F = forces[IND];
Expand All @@ -377,14 +358,20 @@ __global__ void ps_forces(c_number4 *poss, GPU_quat *orientations, c_number4 *fo
CUDA_FS_bond_list bonds;
CUDAStressTensor p_st;

int p_crossover = crossovers[IND];

for(int j = 0; j < num_neighs; j++) {
int q_index = NEXT_NEIGHBOUR(IND, j, matrix_neighs);

if(q_index != IND) {
c_number4 qpos = poss[q_index];
int q_btype = get_particle_btype(qpos);

_WCA(ppos, qpos, MD_WCA_sigma_unbonded[0], MD_sqr_rep_rcut_unbonded[0], F, p_st, box);
bool involve_crossover = p_crossover || crossovers[q_index];
c_number WCA_sigma = (involve_crossover) ? MD_WCA_sigma_crossover[0] : MD_WCA_sigma_unbonded[0];
c_number sqr_rep_rcut = (involve_crossover) ? MD_sqr_rep_rcut_crossover[0] : MD_sqr_rep_rcut_unbonded[0];

_WCA(ppos, qpos, WCA_sigma, sqr_rep_rcut, F, p_st, box);

if(_sticky_interaction(p_btype, q_btype)) {
int eps_idx = p_btype + MD_interaction_matrix_size[0] * q_btype;
Expand Down Expand Up @@ -435,6 +422,9 @@ void CUDACGNucleicAcidsInteraction::cuda_init(int N) {
int tmp_N_strands;
CGNucleicAcidsInteraction::read_topology(&tmp_N_strands, particles);

CUDA_SAFE_CALL(GpuUtils::LR_cudaMalloc(&_d_crossovers, N * sizeof(int)));
CUDA_SAFE_CALL(cudaMemcpy(_d_crossovers, crossovers.data(), N * sizeof(int), cudaMemcpyHostToDevice));

CUDA_SAFE_CALL(GpuUtils::LR_cudaMalloc(&_d_three_body_forces, N * sizeof(c_number4)));
CUDA_SAFE_CALL(GpuUtils::LR_cudaMalloc(&_d_three_body_torques, N * sizeof(c_number4)));

Expand Down Expand Up @@ -492,7 +482,7 @@ void CUDACGNucleicAcidsInteraction::compute_forces(CUDABaseList *lists, c_number

ps_forces
<<<_launch_cfg.blocks, _launch_cfg.threads_per_block>>>
(d_poss, d_orientations, d_forces, d_torques, _d_three_body_forces, _d_three_body_torques, lists->d_matrix_neighs, lists->d_number_neighs, _update_st, _tex_eps, _d_st, d_box);
(d_poss, d_orientations, d_forces, d_torques, _d_three_body_forces, _d_three_body_torques, lists->d_matrix_neighs, lists->d_number_neighs, _d_crossovers, _update_st, _tex_eps, _d_st, d_box);
CUT_CHECK_ERROR("forces_second_step CGNucleicAcids simple_lists error");

// add the three body contributions to the two-body forces and torques
Expand Down
1 change: 1 addition & 0 deletions contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ class CUDACGNucleicAcidsInteraction: public CUDABaseInteraction, public CGNuclei
c_number4 *_d_three_body_forces = nullptr;
c_number4 *_d_three_body_torques = nullptr;
float *_d_3b_epsilon = nullptr;
int *_d_crossovers = nullptr;

cudaTextureObject_t _tex_eps = 0;

Expand Down

0 comments on commit 8ad9ce2

Please sign in to comment.