Merge branch 'master' of https://github.com/lorenzo-rovigatti/oxDNA

analysis/src/oxDNA_analysis_tools/UTILS/boilerplate.py

@@ -1,7 +1,7 @@
 import oxpy
 from oxpy import InputFile
 from os.path import join, abspath, dirname, exists, basename
-from os import mkdir, getpid
+from os import mkdir, getpid,chdir, getcwd
 from multiprocessing import Process
 from shutil import copyfile, rmtree
 from oxDNA_analysis_tools.UTILS.oxview import from_path, oxdna_conf
@@ -14,7 +14,8 @@
 import ipywidgets as widgets
 from IPython.display import display, IFrame
 import numpy as np
-
+from functools import wraps
+import subprocess
 
 default_input_file = {
     "T" :"20C",
@@ -84,7 +85,7 @@ def setup_simulation(top_path:str, dat_path:str, out_dir:str, parameters:dict[st
     out_top = join(out_dir,basename(top_path))
     if(not exists(out_top)):
         copyfile(top_path, out_top)
-    input_file["topology"] = out_top
+    input_file["topology"] = basename(top_path)
 
     #copy the conf over if not present
     out_dat = join(out_dir, basename(dat_path))
@@ -94,12 +95,12 @@ def setup_simulation(top_path:str, dat_path:str, out_dir:str, parameters:dict[st
         copyfile(dat_path, out_dat)
         copyfile(dat_path, last_conf_path)
 
-    input_file["conf_file"] = out_dat
+    input_file["conf_file"] = basename(dat_path)
 
     #set the outputs to desired location
-    input_file["trajectory_file"] = join(out_dir, "trajectory.dat")
-    input_file["energy_file"] = join(out_dir, "energy.dat")
-    input_file["lastconf_file"] =last_conf_path
+    input_file["trajectory_file"] = "trajectory.dat"
+    input_file["energy_file"] = "energy.dat"
+    input_file["lastconf_file"] ="last_conf.dat"
 
     # if we have defined forces we need to 
     # 1) build the file
@@ -109,7 +110,7 @@ def setup_simulation(top_path:str, dat_path:str, out_dir:str, parameters:dict[st
         dump_json(force_dict, force_path)
         input_file["external_forces_as_JSON"] = "true"
         input_file["external_forces"] = "1"
-        input_file["external_forces_file"] = force_path
+        input_file["external_forces_file"] = "forces.json"
 
 
     # set all base simulation parameters, overriding previous
@@ -118,15 +119,38 @@ def setup_simulation(top_path:str, dat_path:str, out_dir:str, parameters:dict[st
 
 
     #prevents output
-    input_file["log_file"] = join(out_dir,"logfile")
+    input_file["log_file"] = "logfile"
 
     input_file_path = join(out_dir,"input")
     with open(input_file_path,"w") as file:
         file.write(str(input_file))
-
-
     return input_file_path
-
+
+# our typical run function for a provided dictionary input file 
+def _prun(input_file_path:str):
+    with oxpy.Context():
+        # change to the directory of the input file path
+        chdir(dirname(abspath(input_file_path)))
+        input_file = InputFile()
+        input_file.init_from_filename(basename(input_file_path))
+        manager = oxpy.OxpyManager(input_file)
+        manager.run_complete() #run complete run's it till the number steps specified are reached 
+
+
+def path_decorator(func):
+    @wraps(func)
+    def wrapper(self, *args, **kwargs):
+        # Save the current path
+        old_path = getcwd()
+        # Change the path to self.__out_dir
+        chdir(self.out_dir)
+        # Call the original method
+        result = func(self, *args, **kwargs)
+        # Restore the original path
+        chdir(old_path)
+        return result
+    return wrapper
+
 
 class Simulation:
     def __init__(self, input_file_path:str):
@@ -138,8 +162,10 @@ def __init__(self, input_file_path:str):
         self.input_file = InputFile()
         self.input_file.init_from_filename(input_file_path)
         self.p = None # the process refference
-        self.out_dir = dirname(input_file_path)
+        self.out_dir = dirname(abspath(input_file_path))
+        self.__input_file_path = input_file_path
 
+    @path_decorator
     def get_init_conf(self):
         """
             returns the initial configuration of the simulation as a rye reader object
@@ -148,6 +174,7 @@ def get_init_conf(self):
                           abspath(self.input_file["conf_file"]))
         return (ti, di), get_confs(ti, di, 0, 1)[0]
 
+    @path_decorator
     def get_last_conf(self):
         """
             returns the last configuration of the simulation as a rye reader object
@@ -171,6 +198,7 @@ def view_last(self):
         (ti,di), conf = self.get_last_conf()
         oxdna_conf(ti, conf)
 
+    @path_decorator
     def get_conf_count(self):
         """
             returns the number of configurations in the trajectory
@@ -179,6 +207,7 @@ def get_conf_count(self):
                          abspath(self.input_file["trajectory_file"]))
         return len(di.idxs)
 
+    @path_decorator
     def get_conf(self, id:int):
         """
             returns the configuration at the given index in the trajectory
@@ -199,6 +228,7 @@ def view_conf(self, id:int):
         (ti,di), conf = self.get_conf(id)
         oxdna_conf(ti, conf)
 
+
     def view_traj(self,  init = 0, op=None):
         """
             opens the trajectory in an embeded oxDNA viewer window
@@ -211,10 +241,11 @@ def view_traj(self,  init = 0, op=None):
 
         slider = widgets.IntSlider(
             min = 0,
-            max = len(di.idxs),
+            max = len(di.idxs)-1,
             step=1,
             description="Select:",
-            value=init
+            value=init, 
+            continuous_update=False,
         )
 
         output = widgets.Output()
@@ -252,22 +283,34 @@ def run(self):
         """
             runs the simulation
         """
-        # our typical run function for a provided dictionary input file 
-        def _prun(input_file:InputFile):
-            with oxpy.Context():
-                manager = oxpy.OxpyManager(input_file)
-                manager.run_complete() #run complete run's it till the number steps specified are reached 
-
+
         # if we have any reference to a simulation
         if(self.p):
             self.terminate()
 
         #spawn the process
-        self.p = Process(target=_prun, args = (self.input_file,))
+        self.p = Process(target=_prun, args = (self.__input_file_path,))
 
         self.p.start()
         return self
 
+    # @path_decorator
+    # TODO: possibly add this to the api
+    # def align(self, p=2):
+    #     """
+    #         aligns the trajectory to the last configuration
+    #     """
+    #     with oxpy.Context():
+    #         input_file = InputFile()
+    #         input_file.init_from_filename(basename(self.__input_file_path))
+    #         print(input_file["trajectory_file"])
+
+    #         # run shell command oat align 
+    #         subprocess.run(["oat", "align", "-p", str(p), input_file["trajectory_file"], "aligned.dat"])
+    #         with open(basename(self.__input_file_path),"w+") as file:
+    #             file.write("aligned = aligned.dat")
+
+    @path_decorator
     def plot_energy(self):
         """
             plots the energy graph of the running simulation
@@ -293,4 +336,4 @@ def terminate(self):
         """
         if(self.p):
             self.p.terminate()
-            print("you evil...")
+            print("you evil...")

contrib/rovigatti/src/Interactions/CUDADetailedPatchySwapInteraction.cu

@@ -494,19 +494,6 @@ void CUDADetailedPatchySwapInteraction::cuda_init(int N) {
 	CUDA_SAFE_CALL(cudaMemcpy(_d_patchy_eps, h_patchy_eps.data(), _patchy_eps.size() * sizeof(float), cudaMemcpyHostToDevice));
 	GpuUtils::init_texture_object(&_tex_patchy_eps, cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat), _d_patchy_eps, _patchy_eps.size());
 
-	int N_base_patches = MAX_PATCHES * N_species;
-	CUDA_SAFE_CALL(GpuUtils::LR_cudaMalloc(&_d_base_patches, N_base_patches * sizeof(float4)));
-	std::vector<float4> h_base_patches(N_base_patches, make_float4(0., 0., 0., 0.));
-	for(uint ns = 0; ns < N_species; ns++) {
-		for(uint np = 0; np < _base_patches[ns].size(); np++) {
-			float4 bp_f4 = make_float4(_base_patches[ns][np].x, _base_patches[ns][np].y, _base_patches[ns][np].z, 0.);
-			h_base_patches[ns * MAX_PATCHES + np] = bp_f4;
-		}
-	}
-
-	CUDA_SAFE_CALL(cudaMemcpy(_d_base_patches, h_base_patches.data(), N_base_patches * sizeof(float4), cudaMemcpyHostToDevice));
-	GpuUtils::init_texture_object(&_tex_base_patches, cudaCreateChannelDesc(32, 32, 32, 32, cudaChannelFormatKindFloat), _d_base_patches, h_base_patches.size());
-
 	for(int i = 0; i < N_species; i++) {
 		int n_patches = _N_patches[i];
 
@@ -522,6 +509,19 @@ void CUDADetailedPatchySwapInteraction::cuda_init(int N) {
 		// fourth argument is the offset
 		CUDA_SAFE_CALL(cudaMemcpyToSymbol(MD_patch_types, patch_types, sizeof(int) * n_patches, i * sizeof(int) * MAX_PATCHES));
 	}
+
+	int N_base_patches = MAX_PATCHES * N_species;
+	CUDA_SAFE_CALL(GpuUtils::LR_cudaMalloc(&_d_base_patches, N_base_patches * sizeof(float4)));
+	std::vector<float4> h_base_patches(N_base_patches, make_float4(0., 0., 0., 0.));
+	for(uint ns = 0; ns < N_species; ns++) {
+		for(uint np = 0; np < _base_patches[ns].size(); np++) {
+			float4 bp_f4 = make_float4(_base_patches[ns][np].x, _base_patches[ns][np].y, _base_patches[ns][np].z, 0.);
+			h_base_patches[ns * MAX_PATCHES + np] = bp_f4;
+		}
+	}
+
+	CUDA_SAFE_CALL(cudaMemcpy(_d_base_patches, h_base_patches.data(), N_base_patches * sizeof(float4), cudaMemcpyHostToDevice));
+	GpuUtils::init_texture_object(&_tex_base_patches, cudaCreateChannelDesc(32, 32, 32, 32, cudaChannelFormatKindFloat), _d_base_patches, h_base_patches.size());
 }
 
 void CUDADetailedPatchySwapInteraction::compute_forces(CUDABaseList *lists, c_number4 *d_poss, GPU_quat *d_orientations, c_number4 *d_forces, c_number4 *d_torques, LR_bonds *d_bonds, CUDABox *d_box) {

contrib/rovigatti/src/Interactions/CUDADetailedPatchySwapInteraction.h

@@ -32,7 +32,7 @@ class CUDADetailedPatchySwapInteraction: public CUDABaseInteraction, public Deta
 
 	llint _step;
 public:
-	static const int MAX_PATCHES = 5;
+	static const int MAX_PATCHES = 20;
 	static const int MAX_SPECIES = 60;
 	static const int MAX_NEIGHS = 20;
 

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/Observables/ObservableOutput.cpp

@@ -119,7 +119,7 @@ void ObservableOutput::_open_output() {
 }
 
 void ObservableOutput::init() {
-	if(!_update_name_with_time) {
+	if(!_update_name_with_time && !_only_last) {
 		_open_output();
 	}
 

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);
 		}