diff --git a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp index 81929584..4e0dd422 100644 --- a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp +++ b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp @@ -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; @@ -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); @@ -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.; @@ -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.; @@ -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; @@ -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; @@ -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; @@ -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)) { @@ -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; @@ -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 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 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(); } @@ -670,6 +660,33 @@ void CGNucleicAcidsInteraction::read_topology(int *N_strands, std::vector(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()) { @@ -767,7 +784,7 @@ void CGNucleicAcidsInteraction::read_topology(int *N_strands, std::vector crossovers; std::map > _bonds; int _N_chains = -1; diff --git a/contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.cu b/contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.cu index 1e520833..11df1dec 100644 --- a/contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.cu +++ b/contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.cu @@ -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]; @@ -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; @@ -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; @@ -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(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(p_st, b2.r, force); -// F += force; -// LR_atomicAddXYZ(forces + b2.q, -force); -// } } } } @@ -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); @@ -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]; @@ -377,6 +358,8 @@ __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); @@ -384,7 +367,11 @@ __global__ void ps_forces(c_number4 *poss, GPU_quat *orientations, c_number4 *fo 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; @@ -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))); @@ -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 diff --git a/contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.h b/contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.h index c68dce22..026c34fe 100644 --- a/contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.h +++ b/contrib/tostiguerra/src/CUDACGNucleicAcidsInteraction.h @@ -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;