From e1eee7c009ed9d19ded16ca84d97b2734e309dbf Mon Sep 17 00:00:00 2001
From: tgfrancesco <57915773+tgfrancesco@users.noreply.github.com>
Date: Wed, 13 Mar 2024 15:35:56 +0100
Subject: [PATCH] fixed rcut for unbonded pair

---
 .../tostiguerra/src/CGNucleicAcidsInteraction.cpp | 15 ++++++++++-----
 .../tostiguerra/src/CGNucleicAcidsInteraction.h   |  2 +-
 2 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp
index 2be3ab33f..fcb4b0a14 100644
--- a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp
+++ b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp
@@ -80,8 +80,9 @@ void CGNucleicAcidsInteraction::get_settings(input_file &inp) {
 void CGNucleicAcidsInteraction::init() {
 	_sqr_rfene = SQR(_rfene);
 	_PS_sqr_rep_rcut = pow(2. * _WCA_sigma, 2. / _PS_n);
-	// _WCA_sigma_unbonded = _WCA_sigma * (6.0 / _bead_size - _3b_sigma) / 2.0; // disabled for now
-	_WCA_sigma_unbonded = _WCA_sigma;
+	_WCA_sigma_unbonded = _WCA_sigma * (6.0 / _bead_size - _3b_sigma) / 2.0; // not disabled
+	_PS_sqr_rep_rcut_unbonded = pow(2. * _WCA_sigma_unbonded, 2. / _PS_n); 
+	// _WCA_sigma_unbonded = _WCA_sigma; // disabled
 
 	OX_LOG(Logger::LOG_INFO, "CGNA: WCA sigma = %lf, WCA sigma unbonded = %lf", _WCA_sigma, _WCA_sigma_unbonded);
 
@@ -215,7 +216,11 @@ 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();
-	if(sqr_r > _PS_sqr_rep_rcut) {
+	
+	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;
+
+	if(sqr_r > sqr_rcut) {
 		return (number) 0.;
 	}
 
@@ -223,10 +228,10 @@ 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;
 
-	number sigma = p->is_bonded(q) ? _WCA_sigma : _WCA_sigma_unbonded;
+	
 
 	// cut-off for all the repulsive interactions
-	if(sqr_r < _PS_sqr_rep_rcut) {
+	if(sqr_r < sqr_rcut) {
 		number part = 1.;
 		number ir2_scaled = SQR(sigma) / sqr_r;
 		for(int i = 0; i < _PS_n / 2; i++) {
diff --git a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.h b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.h
index 091d0e3b4..4fbe99372 100644
--- a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.h
+++ b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.h
@@ -47,7 +47,7 @@ class CGNucleicAcidsInteraction: public BaseInteraction {
 	number _rfene = 1.5;
 	number _sqr_rfene;
 	number _WCA_sigma = 1.0, _WCA_sigma_unbonded;
-	number _PS_sqr_rep_rcut;
+	number _PS_sqr_rep_rcut, _PS_sqr_rep_rcut_unbonded;
 	number _tC = 37.0;
 	number dS_mod = 1.0;
 	number alpha_mod = 1.0;