some changes

src/capabilities/simplemd.cpp

@@ -168,6 +168,7 @@
     m_T0 = Json2KeyWord<double>(m_defaults, "T");
     m_rmrottrans = Json2KeyWord<int>(m_defaults, "rmrottrans");
     m_nocenter = Json2KeyWord<bool>(m_defaults, "nocenter");
+    m_COM = Json2KeyWord<bool>(m_defaults, "COM");
     m_dump = Json2KeyWord<int>(m_defaults, "dump");
     m_print = Json2KeyWord<int>(m_defaults, "print");
     m_max_top_diff = Json2KeyWord<int>(m_defaults, "MaxTopoDiff");
@@ -339,7 +340,7 @@
     m_molecule.setCharge(0);
     if (!m_nocenter) {
         std::cout << "Move stucture to the origin ... " << std::endl;
-        m_molecule.setGeometry(GeometryTools::TranslateGeometry(m_molecule.getGeometry(), GeometryTools::Centroid(m_molecule.getGeometry()), Position{ 0, 0, 0 }));
+        m_molecule.Center(m_COM);
     } else
         std::cout << "Move stucture NOT to the origin ... " << std::endl;
 
@@ -508,7 +509,7 @@
 
 void SimpleMD::InitConstrainedBonds()
 {
 
     if (m_rattle) {
         auto m = m_molecule.DistanceMatrix();
         m_topo_initial = m.second;
@@ -714,6 +715,7 @@
     restart["gradient"] = Tools::DoubleVector2String(m_gradient);
     restart["rmrottrans"] = m_rmrottrans;
     restart["nocenter"] = m_nocenter;
+    restart["COM"] = m_COM;
     restart["average_T"] = m_aver_Temp;
     restart["average_Epot"] = m_aver_Epot;
     restart["average_Ekin"] = m_aver_Ekin;
@@ -819,6 +821,10 @@
         m_nocenter = state["nocenter"];
     } catch (json::type_error& e) {
     }
+    try {
+        m_COM = state["COM"];
+    } catch (json::type_error& e) {
+    }
     try {
         m_T0 = state["T"];
     } catch (json::type_error& e) {
@@ -1109,6 +1115,12 @@
             }
         }
 
+/////////// Dipole calc
+
+
+
+
+/////////// Dipole calc
         if (m_step % m_dump == 0) {
             bool write = WriteGeometry();
             if (write) {
@@ -1518,62 +1530,62 @@
     EKin();
 }
 
 void SimpleMD::ApplyRMSDMTD(double* grad)
 {
     std::chrono::time_point<std::chrono::system_clock> m_start, m_end;
     m_start = std::chrono::system_clock::now();
     m_colvar_incr = 0;
 
     Geometry current_geometry = m_rmsd_mtd_molecule.getGeometry();
     for (int i = 0; i < m_rmsd_indicies.size(); ++i) {
         current_geometry(i, 0) = m_current_geometry[3 * m_rmsd_indicies[i] + 0];
         current_geometry(i, 1) = m_current_geometry[3 * m_rmsd_indicies[i] + 1];
         current_geometry(i, 2) = m_current_geometry[3 * m_rmsd_indicies[i] + 2];
     }
 
     double current_bias = 0;
     double rmsd_reference = 0;
 
     if (m_bias_structure_count == 0) {
         m_bias_threads[0]->addGeometry(current_geometry, 0, m_currentStep, 0);
         m_bias_structure_count++;
         m_rmsd_mtd_molecule.writeXYZFile(Basename() + ".mtd.xyz");
         std::ofstream colvarfile;
         colvarfile.open("COLVAR");
         colvarfile.close();
     }
     if (m_threads == 1 || m_bias_structure_count == 1) {
         for (int i = 0; i < m_bias_threads.size(); ++i) {
             m_bias_threads[i]->setCurrentGeometry(current_geometry, m_currentStep);
             m_bias_threads[i]->start();
             current_bias += m_bias_threads[i]->BiasEnergy();
             for (int j = 0; j < m_rmsd_indicies.size(); ++j) {
                 grad[3 * m_rmsd_indicies[j] + 0] += m_bias_threads[i]->Gradient()(j, 0);
                 grad[3 * m_rmsd_indicies[j] + 1] += m_bias_threads[i]->Gradient()(j, 1);
                 grad[3 * m_rmsd_indicies[j] + 2] += m_bias_threads[i]->Gradient()(j, 2);
             }
             m_colvar_incr += m_bias_threads[i]->Counter();
             m_loop_time += m_bias_threads[i]->Time();
         }
     } else {
         if (m_bias_structure_count < m_threads) {
             for (int i = 0; i < m_bias_structure_count; ++i) {
                 m_bias_threads[i]->setCurrentGeometry(current_geometry, m_currentStep);
             }
         } else {
             for (int i = 0; i < m_bias_threads.size(); ++i) {
                 m_bias_threads[i]->setCurrentGeometry(current_geometry, m_currentStep);
             }
         }
 
         m_bias_pool->setActiveThreadCount(m_threads);
         m_bias_pool->StaticPool();
         m_bias_pool->StartAndWait();
         m_bias_pool->setWakeUp(m_bias_pool->WakeUp() / 2);
 
         for (int i = 0; i < m_bias_threads.size(); ++i) {
             if (m_bias_threads[i]->Return() == 1) {
 
                 current_bias += m_bias_threads[i]->BiasEnergy();
                 for (int j = 0; j < m_rmsd_indicies.size(); ++j) {
                     grad[3 * m_rmsd_indicies[j] + 0] += m_bias_threads[i]->Gradient()(j, 0);

src/capabilities/simplemd.h

@@ -54,7 +54,7 @@
     BiasThread(const Molecule& reference, const json& rmsdconfig);
     ~BiasThread();
 
     virtual int execute() override;
     inline void addGeometry(const Geometry& geometry, double rmsd, double time, int index)
     {
         BiasStructure str;
@@ -130,6 +130,7 @@
     { "Spin", 0 },
     { "rmrottrans", 0 },
     { "nocenter", false },
+    { "COM", false },
     { "dump", 50 },
     { "print", 1000 },
     { "unique", false },
@@ -307,6 +308,7 @@
     bool m_initialised = false, m_restart = false, m_writeUnique = true, m_opt = false, m_rescue = false, m_writeXYZ = true, m_writeinit = false, m_norestart = false;
     int m_rmrottrans = 0, m_rattle_maxiter = 100;
     bool m_nocenter = false;
+    bool m_COM = false;
     bool m_wall_render = false;
     EnergyCalculator* m_interface;
     RMSDTraj* m_unqiue;

src/main.cpp

@@ -893,7 +893,7 @@
             int count = 1;
             double sum = 0, sum_mass = 0, sqrt_sum = 0, sqrt_sum_mass = 0, hmass = 1;
             for (std::size_t i = 2; i < argc; ++i) {
 
                 if (strcmp(argv[i], "-hmass") == 0) {
                     if (i + 1 < argc)
                         hmass = std::stoi(argv[i + 1]);
@@ -953,9 +953,11 @@
             while (!file.AtEnd()) { // calculation and output dipole moment
                 mol = file.Next(); // load Molecule
                 mol.Center(false);
+
                 EnergyCalculator interface("gfn2", blob); // set method to gfn2-xtb and give
                 interface.setMolecule(mol); // set molecule for calc
                 interface.CalculateEnergy(false, true); // calc energy and charges and dipole moment
+
                 mol.setPartialCharges(interface.Charges()); // calc Partial Charges and give it to mol
                 auto charges = interface.Charges(); // dec and init charges
                 mol.setDipole(interface.Dipole() * au);
@@ -1040,7 +1042,7 @@
                       << std::endl;
 
         } else if (strcmp(argv[1], "-stride") == 0) {
 
             if (argc < 4) {
                 std::cerr << "Please use curcuma to center a structure as follows:\ncurcuma -stride trjectory.xyz 100" << std::endl;
                 return 0;