Make the external_force observable work

docs/source/observables.md

@@ -171,6 +171,13 @@ Print the energy associated to all (or a subset of) the external forces acting o
 * `type = force_energy`: the observable type.
 * `[print_group = <string>]`: limit the energy computation to the forces belonging to a specific group of forces. This can be set by adding a `group_name` option to the [desired external forces](forces.md#common-options). If not set, all external forces will be considered.
 
+## External force acting on particle(s)
+
+Print the force vector acting on all (or a subset of all) particles due to external forces. This observable supports the `update_every` option.
+
+* `type = external_force`: the observable type.
+* `particles`: list of comma-separated particle indexes whose force vectors should be printed.
+
 ## Configuration
 
 Print an [oxDNA configuration](configurations.md#configuration-file).

src/Backends/FIREBackend.cpp

@@ -162,7 +162,7 @@ void FIREBackend::_first_step() {
 			p->torque = LR_vector((number) 0, (number) 0, (number) 0);
 
 		}
-		p->set_initial_forces(current_step(), _box);
+		p->set_initial_forces(current_step(), _box.get());
 
 		_lists->single_update(p);
 	}

src/Backends/MD_CPUBackend.cpp

@@ -128,7 +128,7 @@ void MD_CPUBackend::_first_step() {
 			p->set_positions();
 		}
 
-		p->set_initial_forces(current_step(), _box);
+		p->set_initial_forces(current_step(), _box.get());
 
 		_lists->single_update(p);
 	}

src/Backends/MinBackend.cpp

@@ -115,7 +115,7 @@ void MinBackend::sim_step() {
 	_mytimer->resume();
 
 	for(auto p: _particles) {
-		p->set_initial_forces(current_step(), _box);
+		p->set_initial_forces(current_step(), _box.get());
 	}
 
 	_timer_lists->resume();

src/Observables/ExternalForce.cpp

@@ -18,30 +18,43 @@ ExternalForce::~ExternalForce() {
 void ExternalForce::get_settings(input_file &my_inp, input_file &sim_inp) {
 	BaseObservable::get_settings(my_inp, sim_inp);
 
-	std::string ids;
-	if(getInputString(&my_inp, "ids", ids, 0) == KEY_FOUND) {
-		_ids = Utils::split(ids, ',');
-	}
+	getInputString(&my_inp, "particles", _particle_ids, 0);
 }
 
-void ExternalForce::update_data(llint curr_step) {
-	static bool initialised = false;
-	if(!initialised) {
-		if(_ids.size() == 0) {
-			_force_averages.resize(_config_info->forces.size());
-			_ext_forces = _config_info->forces;
-		}
-		else {
-			_force_averages.resize(_ids.size());
-			for(auto id : _ids) {
-				_ext_forces.push_back(_config_info->get_force_by_id(id));
+void ExternalForce::init() {
+	std::vector<int> ids = Utils::get_particles_from_string(_config_info->particles(), _particle_ids, "ExternalForce observable");
+
+	if(ids.size() == 0) {
+		_particles = _config_info->particles();
+	}
+	else {
+		for(auto p : _config_info->particles()) {
+			if(std::find(ids.begin(), ids.end(), p->index) != std::end(ids)) {
+				_particles.push_back(p);
 			}
-			_ext_forces.resize(_ids.size());
 		}
-		initialised = true;
 	}
+	_force_averages.resize(_particles.size());
+}
+
+void ExternalForce::update_data(llint curr_step) {
+	for(uint i = 0; i < _particles.size(); i++) {
+		auto p = _particles[i];
+		p->set_initial_forces(curr_step, _config_info->box);
+		_force_averages[i] += p->force;
+	}
+	_times_updated++;
 }
 
 std::string ExternalForce::get_output_string(llint curr_step) {
-	return Utils::sformat("TO BE IMPLEMENTED");
+	std::string output = Utils::sformat("# t = %lld\n", curr_step);
+
+	for(auto &force : _force_averages) {
+		force /= _times_updated;
+		output += Utils::sformat("%lf %lf %lf\n", force.x, force.y, force.z);
+		force = LR_vector();
+	}
+	_times_updated = 0;
+
+	return output;
 }

src/Observables/ExternalForce.h

@@ -16,14 +16,16 @@
 
 class ExternalForce: public BaseObservable {
 protected:
-	std::vector<std::string> _ids;
-	std::vector<std::shared_ptr<BaseForce>> _ext_forces;
+	std::string _particle_ids;
+	std::vector<BaseParticle *> _particles;
 	std::vector<LR_vector> _force_averages;
 public:
 	ExternalForce();
 	virtual ~ExternalForce();
 
-	virtual void get_settings(input_file &my_inp, input_file &sim_inp);
+	void get_settings(input_file &my_inp, input_file &sim_inp) override;
+	void init() override;
+
 	void update_data(llint curr_step) override;
 
 	std::string get_output_string(llint curr_step);

src/Particles/BaseParticle.h

@@ -73,7 +73,7 @@ class BaseParticle {
 	 */
 	bool add_ext_force(BaseForce *f);
 
-	inline void set_initial_forces(llint step, const std::shared_ptr<BaseBox> &box) {
+	inline void set_initial_forces(llint step, BaseBox *box) {
 		if(is_rigid_body()) {
 			torque = LR_vector((number) 0.f, (number) 0.f, (number) 0.f);
 		}