Take into account IClamp for the lfp calculation

Author: jorblancoa

coreneuron/io/reports/report_event.cpp

@@ -80,7 +80,7 @@ void ReportEvent::lfp_calc(NrnThread* nt) {
     if (step > 0 && (static_cast<int>(step) % reporting_period) == 0) {
         auto* mapinfo = static_cast<NrnThreadMappingInfo*>(nt->mapping);
         double* fast_imem_rhs = nt->nrn_fast_imem->nrn_sav_rhs;
-
+        auto& summation_report = nt->summation_report_handler_->summation_reports_[report_path];
         for (const auto& kv: vars_to_report) {
             int gid = kv.first;
             const auto& to_report = kv.second;
@@ -94,7 +94,13 @@ void ReportEvent::lfp_calc(NrnThread* nt) {
                     if(std::isnan(factor)) {
                         factor = 0.0;
                     }
-                    sum += fast_imem_rhs[segment_id] * factor;
+                    double iclamp = 0.0;
+                    for (const auto& value: summation_report.currents_[segment_id]) {
+                        double current_value = *value.first;
+                        int scale = value.second;
+                        iclamp += current_value * scale;
+                    }
+                    sum += (fast_imem_rhs[segment_id] + iclamp) * factor;
                     count++;
             }
             *(to_report.front().var_value) = sum;
@@ -107,8 +113,10 @@ void ReportEvent::deliver(double t, NetCvode* nc, NrnThread* nt) {
 /* reportinglib is not thread safe */
 #pragma omp critical
     {
-        summation_alu(nt);
-        if (report_type == ReportType::LFPReport) {
+        if (report_type == ReportType::SummationReport) {
+            summation_alu(nt);
+        }
+        else if (report_type == ReportType::LFPReport) {
             lfp_calc(nt);
         }
         // each thread needs to know its own step

coreneuron/io/reports/report_handler.cpp

@@ -62,7 +62,7 @@ void ReportHandler::create_report(double dt, double tstop, double delay) {
             case LFPReport:
                 // 1 lfp value per gid
                 mapinfo->_lfp.resize(nt.ncell);
-                vars_to_report = get_lfp_vars_to_report(nt, m_report_config, mapinfo->_lfp.data());
+                vars_to_report = get_lfp_vars_to_report(nt, m_report_config, mapinfo->_lfp.data(), nodes_to_gid);
                 is_soma_target = m_report_config.section_type == SectionType::Soma ||
                                  m_report_config.section_type == SectionType::Cell;
                 register_section_report(nt, m_report_config, vars_to_report, is_soma_target);
@@ -353,14 +353,25 @@ VarsToReport ReportHandler::get_synapse_vars_to_report(
 
 VarsToReport ReportHandler::get_lfp_vars_to_report(const NrnThread& nt,
                                                    ReportConfiguration& report,
-                                                   double* report_variable) const {
+                                                   double* report_variable,
+                                                   const std::vector<int>& nodes_to_gids) const {
+    auto& summation_report = nt.summation_report_handler_->summation_reports_[report.output_path];
     VarsToReport vars_to_report;
     for (int i = 0; i < nt.ncell; i++) {
         int gid = nt.presyns[i].gid_;
         if (report.target.find(gid) == report.target.end()) {
             continue;
         }
-
+        // IClamp is needed for the LFP calculation
+        auto mech_id = nrn_get_mechtype("IClamp");
+        Memb_list* ml = nt._ml_list[mech_id];
+        for (int j = 0; j < ml->nodecount; j++) {
+            auto segment_id = ml->nodeindices[j];
+            if ((nodes_to_gids[segment_id] == gid)) {
+                double* var_value = get_var_location_from_var_name(mech_id, "i", ml, j);
+                summation_report.currents_[segment_id].push_back(std::make_pair(var_value, -1));
+            }
+        }
         std::vector<VarWithMapping> to_report;
         double* variable = report_variable + i;
         to_report.push_back(VarWithMapping(i, variable));

coreneuron/io/reports/report_handler.hpp

@@ -46,7 +46,8 @@ class ReportHandler {
                                             const std::vector<int>& nodes_to_gids) const;
     VarsToReport get_lfp_vars_to_report(const NrnThread& nt,
                                         ReportConfiguration& report,
-                                        double* report_variable) const;
+                                        double* report_variable,
+                                        const std::vector<int>& nodes_to_gids) const;
     std::vector<int> map_gids(const NrnThread& nt) const;
 #endif  // defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS)
   protected: