[e9cfc4] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2012 University of Bonn. All rights reserved.
|
---|
[5aaa43] | 5 | * Copyright (C) 2013 Frederik Heber. All rights reserved.
|
---|
[e9cfc4] | 6 | *
|
---|
| 7 | *
|
---|
| 8 | * This file is part of MoleCuilder.
|
---|
| 9 | *
|
---|
| 10 | * MoleCuilder is free software: you can redistribute it and/or modify
|
---|
| 11 | * it under the terms of the GNU General Public License as published by
|
---|
| 12 | * the Free Software Foundation, either version 2 of the License, or
|
---|
| 13 | * (at your option) any later version.
|
---|
| 14 | *
|
---|
| 15 | * MoleCuilder is distributed in the hope that it will be useful,
|
---|
| 16 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 17 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 18 | * GNU General Public License for more details.
|
---|
| 19 | *
|
---|
| 20 | * You should have received a copy of the GNU General Public License
|
---|
| 21 | * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
|
---|
| 22 | */
|
---|
| 23 |
|
---|
| 24 | /*
|
---|
| 25 | * InterfaceVMGJob.cpp
|
---|
| 26 | *
|
---|
| 27 | * Created on: 10.06.2012
|
---|
| 28 | * Author: Frederik Heber
|
---|
| 29 | */
|
---|
| 30 |
|
---|
| 31 | #ifdef HAVE_CONFIG_H
|
---|
| 32 | #include <config.h>
|
---|
| 33 | #endif
|
---|
| 34 |
|
---|
[cd77fc] | 35 | #ifdef HAVE_MPI
|
---|
| 36 | #include "mpi.h"
|
---|
| 37 | #endif
|
---|
| 38 |
|
---|
[e9cfc4] | 39 | #include "base/vector.hpp"
|
---|
[cd77fc] | 40 | #include "base/math.hpp"
|
---|
[17fcbe7] | 41 | #include "comm/comm.hpp"
|
---|
[e9cfc4] | 42 | #include "grid/grid.hpp"
|
---|
| 43 | #include "grid/multigrid.hpp"
|
---|
[a82602] | 44 | #include "units/particle/comm_mpi_particle.hpp"
|
---|
| 45 | #include "units/particle/interpolation.hpp"
|
---|
| 46 | #include "units/particle/linked_cell_list.hpp"
|
---|
[17fcbe7] | 47 | #include "mg.hpp"
|
---|
[e9cfc4] | 48 |
|
---|
| 49 | #include "InterfaceVMGJob.hpp"
|
---|
| 50 |
|
---|
[0990e1] | 51 | #include "CodePatterns/MemDebug.hpp"
|
---|
| 52 |
|
---|
| 53 | #include <cmath>
|
---|
| 54 | #include <iostream>
|
---|
[8f3cdd] | 55 | #include <limits>
|
---|
[0990e1] | 56 |
|
---|
| 57 | #include "CodePatterns/Log.hpp"
|
---|
| 58 |
|
---|
[b123a5] | 59 | #include "Jobs/WindowGrid_converter.hpp"
|
---|
[0990e1] | 60 |
|
---|
[e9cfc4] | 61 | using namespace VMG;
|
---|
| 62 | using VMGInterfaces::InterfaceVMGJob;
|
---|
| 63 |
|
---|
[8f3cdd] | 64 | InterfaceVMGJob::InterfaceVMGJob(const SamplingGrid &_sampled_input,
|
---|
[2bc560] | 65 | VMGData &_returndata,
|
---|
[cd77fc] | 66 | const std::vector< std::vector<double> > &_particle_positions,
|
---|
| 67 | const std::vector< double > &_particle_charges,
|
---|
| 68 | VMG::Boundary boundary,
|
---|
| 69 | int levelMin,
|
---|
| 70 | int levelMax,
|
---|
[8f3cdd] | 71 | const VMG::Vector &_box_begin,
|
---|
| 72 | vmg_float _box_end,
|
---|
[cd77fc] | 73 | const int& near_field_cells,
|
---|
[e2925fd] | 74 | const ImportParticles_t _ImportParticles,
|
---|
[b6b21a] | 75 | const bool _DoPrintDebug,
|
---|
[cd77fc] | 76 | int coarseningSteps,
|
---|
| 77 | double alpha) :
|
---|
| 78 | VMG::Interface(boundary, levelMin, levelMax,
|
---|
[8f3cdd] | 79 | _box_begin, _box_end, coarseningSteps, alpha),
|
---|
[cd77fc] | 80 | spl(near_field_cells, Extent(MaxLevel()).MeshWidth().Max()),
|
---|
| 81 | sampled_input(_sampled_input),
|
---|
[2bc560] | 82 | returndata(_returndata),
|
---|
[b6b21a] | 83 | level(levelMax),
|
---|
[e2925fd] | 84 | ImportParticles(_ImportParticles),
|
---|
[b47f9c] | 85 | DoPrintDebug(_DoPrintDebug),
|
---|
| 86 | OpenBoundaryCondition(boundary[0] == VMG::Open)
|
---|
[cd77fc] | 87 | {
|
---|
[8f3cdd] | 88 | for (size_t i=0;i<3;++i) {
|
---|
| 89 | box_begin[i] = _box_begin[i];
|
---|
| 90 | box_end[i] = _box_end;
|
---|
| 91 | }
|
---|
[cd77fc] | 92 | std::vector< std::vector<double> >::const_iterator positer = _particle_positions.begin();
|
---|
| 93 | std::vector<double>::const_iterator chargeiter = _particle_charges.begin();
|
---|
| 94 | double pos[3];
|
---|
| 95 | for (; positer != _particle_positions.end(); ++positer, ++chargeiter) {
|
---|
| 96 | ASSERT( (*positer).size() == 3,
|
---|
| 97 | "InterfaceVMGJob::InterfaceVMGJob() - particle "
|
---|
| 98 | +toString(distance(_particle_positions.begin(), positer))+" has not exactly 3 coordinates.");
|
---|
| 99 | for (size_t i=0;i<3;++i)
|
---|
| 100 | pos[i] = (*positer)[i];
|
---|
| 101 | particles.push_back(Particle::Particle(pos, *chargeiter));
|
---|
| 102 | }
|
---|
| 103 | }
|
---|
| 104 |
|
---|
[e9cfc4] | 105 | void InterfaceVMGJob::ImportRightHandSide(Multigrid& multigrid)
|
---|
| 106 | {
|
---|
| 107 | Index i;
|
---|
| 108 | Vector pos;
|
---|
[cd77fc] | 109 | // VMG::TempGrid *temp_grid = new VMG::TempGrid(129, 0, 0., 1.);
|
---|
[e9cfc4] | 110 |
|
---|
| 111 | Grid& grid = multigrid(multigrid.MaxLevel());
|
---|
[a82602] | 112 | grid.Clear();
|
---|
[17fcbe7] | 113 | //grid.ClearBoundary(); // we don't have a boundary under periodic boundary conditions
|
---|
[e9cfc4] | 114 |
|
---|
[a1436b] | 115 | // print debugging info on grid size
|
---|
| 116 | LOG(1, "INFO: Mesh has extent " << grid.Extent().MeshWidth() << ".");
|
---|
| 117 | const int gridpoints = pow(2, level);
|
---|
| 118 | LOG(1, "INFO: gridpoints on finest level are " << gridpoints << ".");
|
---|
| 119 | LOG(1, "INFO: "
|
---|
| 120 | << "X in [" << grid.Local().Begin().X() << "," << grid.Local().End().X() << "],"
|
---|
| 121 | << "Y in [" << grid.Local().Begin().Y() << "," << grid.Local().End().Y() << "],"
|
---|
| 122 | << "Z in [" << grid.Local().Begin().Z() << "," << grid.Local().End().Z() << "].");
|
---|
| 123 |
|
---|
[cd77fc] | 124 | /// 1. assign nuclei as smeared-out charges to the grid
|
---|
| 125 |
|
---|
| 126 | /*
|
---|
| 127 | * Charge assignment on the grid
|
---|
| 128 | */
|
---|
[a82602] | 129 | Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
|
---|
[cd77fc] | 130 | Grid& particle_grid = comm.GetParticleGrid();
|
---|
| 131 | particle_grid.Clear();
|
---|
| 132 |
|
---|
[a82602] | 133 | // distribute particles
|
---|
| 134 | particles.clear();
|
---|
| 135 | comm.CommParticles(grid, particles);
|
---|
| 136 |
|
---|
[cd77fc] | 137 | assert(particle_grid.Global().LocalSize().IsComponentwiseGreater(
|
---|
| 138 | VMG::MG::GetFactory().GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS")));
|
---|
| 139 |
|
---|
[e2925fd] | 140 | if (ImportParticles == DoImportParticles) {
|
---|
| 141 | // create smeared-out particle charges on particle_grid via splines
|
---|
| 142 | LOG(1, "INFO: Creating particle grid for " << particles.size() << " particles.");
|
---|
| 143 | for (std::list<Particle::Particle>::iterator iter = particles.begin();
|
---|
| 144 | iter != particles.end(); ++iter) {
|
---|
| 145 | LOG(2, "DEBUG: Current particle is at " << (*iter).Pos()
|
---|
| 146 | << " with charge " << (*iter).Charge() << ".");
|
---|
| 147 | spl.SetSpline(particle_grid, *iter);
|
---|
| 148 | }
|
---|
[4adef7] | 149 | }
|
---|
| 150 |
|
---|
[cd77fc] | 151 | // Communicate charges over halo
|
---|
| 152 | comm.CommFromGhosts(particle_grid);
|
---|
| 153 |
|
---|
[b6b21a] | 154 | if (DoPrintDebug) {
|
---|
| 155 | // print nuclei grid to vtk
|
---|
| 156 | comm.PrintGrid(particle_grid, "Sampled Nuclei Density");
|
---|
| 157 | }
|
---|
[17fcbe7] | 158 |
|
---|
| 159 | // add sampled electron charge density onto grid
|
---|
[b123a5] | 160 | WindowGrid_converter::addWindowOntoGrid(
|
---|
[8f3cdd] | 161 | grid,
|
---|
| 162 | sampled_input,
|
---|
[b47f9c] | 163 | 1.,
|
---|
| 164 | OpenBoundaryCondition);
|
---|
[17fcbe7] | 165 |
|
---|
[b6b21a] | 166 | if (DoPrintDebug) {
|
---|
| 167 | // print electron grid to vtk
|
---|
| 168 | comm.PrintGrid(grid, "Sampled Electron Density");
|
---|
| 169 | }
|
---|
[a1436b] | 170 |
|
---|
[17fcbe7] | 171 | // add particle_grid onto grid
|
---|
[cd77fc] | 172 | for (int i=0; i<grid.Local().Size().X(); ++i)
|
---|
| 173 | for (int j=0; j<grid.Local().Size().Y(); ++j)
|
---|
| 174 | for (int k=0; k<grid.Local().Size().Z(); ++k)
|
---|
| 175 | grid(grid.Local().Begin().X() + i,
|
---|
| 176 | grid.Local().Begin().Y() + j,
|
---|
[a1436b] | 177 | grid.Local().Begin().Z() + k) = 4.0 * VMG::Math::pi * (
|
---|
| 178 | grid(grid.Local().Begin().X() + i,
|
---|
| 179 | grid.Local().Begin().Y() + j,
|
---|
| 180 | grid.Local().Begin().Z() + k) +
|
---|
[cd77fc] | 181 | particle_grid.GetVal(particle_grid.Local().Begin().X() + i,
|
---|
| 182 | particle_grid.Local().Begin().Y() + j,
|
---|
[a1436b] | 183 | particle_grid.Local().Begin().Z() + k));
|
---|
[d12d621] | 184 |
|
---|
[17fcbe7] | 185 | // calculate sum over grid times h^3 as check, should be roughly zero
|
---|
| 186 | const double element_volume = grid.Extent().MeshWidth().Product();
|
---|
[e9cfc4] | 187 | double charge_sum = 0.0;
|
---|
[08cc62] | 188 | for (Grid::iterator grid_iter = grid.Iterators().Local().Begin();
|
---|
| 189 | grid_iter != grid.Iterators().Local().End();
|
---|
| 190 | ++grid_iter)
|
---|
[e9cfc4] | 191 | charge_sum += grid.GetVal(*grid_iter);
|
---|
[17fcbe7] | 192 | charge_sum = element_volume * comm.GlobalSum(charge_sum);
|
---|
[e8e472b] | 193 | comm.PrintOnce(Debug, "Grid charge integral: %e", charge_sum/(4.0 * VMG::Math::pi));
|
---|
[e9cfc4] | 194 |
|
---|
[b6b21a] | 195 | if (DoPrintDebug) {
|
---|
| 196 | // print total grid to vtk
|
---|
| 197 | comm.PrintGrid(grid, "Total Charge Density");
|
---|
| 198 | }
|
---|
[16b4fe] | 199 |
|
---|
[e9cfc4] | 200 | // delete temp_grid;
|
---|
| 201 | }
|
---|
| 202 |
|
---|
| 203 | void InterfaceVMGJob::ExportSolution(Grid& grid)
|
---|
| 204 | {
|
---|
[a82602] | 205 | /// sample the obtained potential to evaluate with the electron charge density
|
---|
| 206 |
|
---|
[e9cfc4] | 207 | // grid now contains the sough-for potential
|
---|
[a82602] | 208 | //Comm& comm = *MG::GetComm();
|
---|
| 209 | Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
|
---|
[e9cfc4] | 210 |
|
---|
| 211 |
|
---|
[b6b21a] | 212 | if (DoPrintDebug) {
|
---|
| 213 | // print output grid to vtk
|
---|
| 214 | comm.PrintGrid(grid, "Potential Solution");
|
---|
| 215 | }
|
---|
[16b4fe] | 216 |
|
---|
[17fcbe7] | 217 | // obtain sampled potential from grid
|
---|
[8f3cdd] | 218 | returndata.sampled_potential.setWindow(
|
---|
| 219 | box_begin,
|
---|
| 220 | box_end
|
---|
| 221 | );
|
---|
[b123a5] | 222 | WindowGrid_converter::addGridOntoWindow(
|
---|
[8f3cdd] | 223 | grid,
|
---|
| 224 | returndata.sampled_potential,
|
---|
[b47f9c] | 225 | +1.,
|
---|
| 226 | OpenBoundaryCondition
|
---|
[8f3cdd] | 227 | );
|
---|
[e9cfc4] | 228 |
|
---|
[17fcbe7] | 229 | // calculate integral over potential as long-range energy contribution
|
---|
| 230 | const double element_volume =
|
---|
| 231 | grid.Extent().MeshWidth().X() * grid.Extent().MeshWidth().Y() * grid.Extent().MeshWidth().Z();
|
---|
[e9cfc4] | 232 | Grid::iterator grid_iter;
|
---|
| 233 | double potential_sum = 0.0;
|
---|
| 234 | for (grid_iter=grid.Iterators().Local().Begin(); grid_iter!=grid.Iterators().Local().End(); ++grid_iter)
|
---|
| 235 | potential_sum += grid.GetVal(*grid_iter);
|
---|
[17fcbe7] | 236 | potential_sum = element_volume * comm.GlobalSum(potential_sum);
|
---|
[e8e472b] | 237 | comm.PrintOnce(Debug, "Grid potential sum: %e", potential_sum);
|
---|
[e9cfc4] | 238 |
|
---|
[a82602] | 239 | {
|
---|
| 240 | Grid::iterator grid_iter = grid.Iterators().Local().Begin();
|
---|
[e8e472b] | 241 | comm.PrintOnce(Debug, "Grid potential at (0,0,0): %e", grid.GetVal(*grid_iter));
|
---|
[a82602] | 242 | }
|
---|
| 243 |
|
---|
| 244 | //Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm()); returndata.e_long = potential_sum;
|
---|
| 245 |
|
---|
| 246 | /// Calculate potential energy of nuclei
|
---|
| 247 |
|
---|
| 248 | vmg_float e = 0.0;
|
---|
| 249 | vmg_float e_long = 0.0;
|
---|
| 250 | vmg_float e_self = 0.0;
|
---|
| 251 | vmg_float e_short_peak = 0.0;
|
---|
| 252 | vmg_float e_short_spline = 0.0;
|
---|
| 253 |
|
---|
| 254 | Factory& factory = MG::GetFactory();
|
---|
| 255 |
|
---|
| 256 | /*
|
---|
| 257 | * Get parameters and arrays
|
---|
| 258 | */
|
---|
| 259 | const vmg_int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
|
---|
| 260 | const vmg_int& interpolation_degree = factory.GetObjectStorageVal<int>("PARTICLE_INTERPOLATION_DEGREE");
|
---|
| 261 |
|
---|
| 262 | Particle::Interpolation ip(interpolation_degree);
|
---|
| 263 |
|
---|
| 264 | const vmg_float r_cut = near_field_cells * grid.Extent().MeshWidth().Max();
|
---|
| 265 |
|
---|
| 266 | /*
|
---|
| 267 | * Copy potential values to a grid with sufficiently large halo size.
|
---|
| 268 | * This may be optimized in future.
|
---|
| 269 | * The parameters of this grid have been set in the import step.
|
---|
| 270 | */
|
---|
| 271 | Grid& particle_grid = comm.GetParticleGrid();
|
---|
| 272 |
|
---|
[b47f9c] | 273 | {
|
---|
| 274 | Index i;
|
---|
| 275 | for (i.X()=0; i.X()<grid.Local().Size().X(); ++i.X())
|
---|
| 276 | for (i.Y()=0; i.Y()<grid.Local().Size().Y(); ++i.Y())
|
---|
| 277 | for (i.Z()=0; i.Z()<grid.Local().Size().Z(); ++i.Z())
|
---|
| 278 | particle_grid(i + particle_grid.Local().Begin()) = grid.GetVal(i + grid.Local().Begin());
|
---|
| 279 | comm.CommToGhosts(particle_grid);
|
---|
| 280 | }
|
---|
[a82602] | 281 |
|
---|
| 282 | /*
|
---|
| 283 | * Compute potentials
|
---|
| 284 | */
|
---|
| 285 | Particle::LinkedCellList lc(particles, near_field_cells, grid);
|
---|
| 286 | Particle::LinkedCellList::iterator p1, p2;
|
---|
| 287 | Grid::iterator iter;
|
---|
| 288 |
|
---|
| 289 | comm.CommLCListToGhosts(lc);
|
---|
| 290 |
|
---|
| 291 | for (int i=lc.Local().Begin().X(); i<lc.Local().End().X(); ++i)
|
---|
| 292 | for (int j=lc.Local().Begin().Y(); j<lc.Local().End().Y(); ++j)
|
---|
| 293 | for (int k=lc.Local().Begin().Z(); k<lc.Local().End().Z(); ++k) {
|
---|
| 294 |
|
---|
| 295 | if (lc(i,j,k).size() > 0)
|
---|
| 296 | ip.ComputeCoefficients(particle_grid, Index(i,j,k) - lc.Local().Begin() + particle_grid.Local().Begin());
|
---|
| 297 |
|
---|
| 298 | for (p1=lc(i,j,k).begin(); p1!=lc(i,j,k).end(); ++p1) {
|
---|
| 299 |
|
---|
| 300 | // Interpolate long-range part of potential and electric field
|
---|
| 301 | ip.Evaluate(**p1);
|
---|
| 302 |
|
---|
| 303 | // Subtract self-induced potential
|
---|
| 304 | (*p1)->Pot() -= (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
|
---|
| 305 |
|
---|
| 306 | e_long += 0.5 * (*p1)->Charge() * ip.EvaluatePotentialLR(**p1);
|
---|
| 307 | e_self += 0.5 * (*p1)->Charge() * (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
|
---|
| 308 |
|
---|
| 309 | for (int dx=-1*near_field_cells; dx<=near_field_cells; ++dx)
|
---|
| 310 | for (int dy=-1*near_field_cells; dy<=near_field_cells; ++dy)
|
---|
| 311 | for (int dz=-1*near_field_cells; dz<=near_field_cells; ++dz) {
|
---|
| 312 |
|
---|
| 313 | for (p2=lc(i+dx,j+dy,k+dz).begin(); p2!=lc(i+dx,j+dy,k+dz).end(); ++p2)
|
---|
| 314 |
|
---|
| 315 | if (*p1 != *p2) {
|
---|
| 316 |
|
---|
| 317 | const Vector dir = (*p1)->Pos() - (*p2)->Pos();
|
---|
| 318 | const vmg_float length = dir.Length();
|
---|
| 319 |
|
---|
| 320 | if (length < r_cut) {
|
---|
| 321 |
|
---|
| 322 | (*p1)->Pot() += (*p2)->Charge() / length * (1.0 + spl.EvaluatePotential(length));
|
---|
| 323 | (*p1)->Field() += (*p2)->Charge() * dir * spl.EvaluateField(length);
|
---|
| 324 |
|
---|
| 325 | e_short_peak += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
|
---|
| 326 | e_short_spline += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
|
---|
| 327 | }
|
---|
| 328 | }
|
---|
| 329 | }
|
---|
| 330 | }
|
---|
| 331 | }
|
---|
| 332 |
|
---|
| 333 | /* Remove average force term */
|
---|
[e4721d3] | 334 | if (!particles.empty()) {
|
---|
| 335 | Vector average_force = 0.0;
|
---|
| 336 | for (std::list<Particle::Particle>::const_iterator iter=particles.begin(); iter!=particles.end(); ++iter)
|
---|
| 337 | average_force += iter->Charge() * iter->Field();
|
---|
| 338 | const vmg_int& npl = MG::GetFactory().GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
|
---|
| 339 | const vmg_int num_particles_global = comm.GlobalSum(npl);
|
---|
| 340 | average_force /= num_particles_global;
|
---|
| 341 | comm.GlobalSumArray(average_force.vec(), 3);
|
---|
| 342 | for (std::list<Particle::Particle>::iterator iter=particles.begin(); iter!=particles.end(); ++iter)
|
---|
| 343 | iter->Field() -= average_force / iter->Charge();
|
---|
| 344 | }
|
---|
[a82602] | 345 |
|
---|
| 346 | comm.CommParticlesBack(particles);
|
---|
| 347 |
|
---|
| 348 | vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");
|
---|
| 349 | const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
|
---|
| 350 | const vmg_float* p = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
|
---|
| 351 | // const vmg_float* f = factory.GetObjectStorageArray<vmg_float>("PARTICLE_FIELD_ARRAY");
|
---|
| 352 |
|
---|
| 353 |
|
---|
| 354 | e_long = comm.GlobalSumRoot(e_long);
|
---|
| 355 | e_short_peak = comm.GlobalSumRoot(e_short_peak);
|
---|
| 356 | e_short_spline = comm.GlobalSumRoot(e_short_spline);
|
---|
| 357 | e_self = comm.GlobalSumRoot(e_self);
|
---|
| 358 |
|
---|
| 359 | for (int j=0; j<num_particles_local; ++j)
|
---|
| 360 | e += 0.5 * p[j] * q[j];
|
---|
| 361 | e = comm.GlobalSumRoot(e);
|
---|
| 362 |
|
---|
[e8e472b] | 363 | comm.PrintOnce(Debug, "E_long: %e", e_long);
|
---|
| 364 | comm.PrintOnce(Debug, "E_short_peak: %e", e_short_peak);
|
---|
| 365 | comm.PrintOnce(Debug, "E_short_spline: %e", e_short_spline);
|
---|
| 366 | comm.PrintOnce(Debug, "E_self: %e", e_self);
|
---|
| 367 | comm.PrintOnce(Debug, "E_total: %e", e);
|
---|
| 368 | comm.PrintOnce(Debug, "E_total*: %e", e_long + e_short_peak + e_short_spline - e_self);
|
---|
[a82602] | 369 |
|
---|
[e2925fd] | 370 | returndata.nuclei_long = e_long;
|
---|
| 371 | returndata.electron_long = e_long;
|
---|
[e9cfc4] | 372 | }
|
---|