Merge branch 'normalization_correction'

multiego.py

@@ -152,6 +152,27 @@ def meGO_parsing():
         type=str,
         help="Symmetry file for the system",
     )
+    optional_args.add_argument(
+        "--f",
+        default=1,
+        type=float,
+        help="partition function normalization",
+    )
+    optional_args.add_argument(
+        "--inter_f",
+        default=1,
+        type=float,
+        help="partition function normalization inter-molecular",
+    )
+    optional_args.add_argument(
+        "--inter_domain_f",
+        default=1,
+        type=float,
+        help="partition function normalization inter_domain",
+    )
+    optional_args.add_argument(
+        "--relative_c12d", default=0.001, type=float, help="Relative deviation from default to set new replulsive c12"
+    )
 
     args, remaining = parser.parse_known_args()
     args.root_dir = os.path.dirname(os.path.abspath(__file__))

test/test_inputs/lyso-bnz_ref/aa_sym_all

@@ -0,0 +1,11 @@
+ASP OD1 OD2
+GLU OE1 OE2
+LEU CD1 CD2
+PHE CD1 CD2
+PHE CE1 CE2
+TYR CD1 CD2
+TYR CE1 CE2
+ARG NH1 NH2
+VAL CG1 CG2
+LYS O1 O2
+BNZ CD1 CD2 CE1 CE2 CZ CG 

tools/cmdata/src/cmdata/cmdata.hpp

@@ -122,6 +122,7 @@ class CMData
     const std::function<ftype_cross_::signature> &f_inter_mol_cross_, double weight
   )
   {
+    semaphore_.acquire();
     const char * atomname;
     int tmp_i = 0;
     std::size_t mol_i = i, mol_j = 0;
@@ -133,7 +134,6 @@ class CMData
     }
     if (mol_i == num_mol_unique_[mol_id_[i]]) mol_i = 0;
     int molb = 0;
-    semaphore_.acquire();
     /* Loop over molecules  */
     for (int j = 0; j < nindex_; j++)
     {
@@ -226,6 +226,7 @@ class CMData
     }
     semaphore_.release();
   }
+    PbcType pbcType_;
 
 public:
   CMData(
@@ -238,7 +239,6 @@ class CMData
       no_pbc_(no_pbc), dt_(dt), t_begin_(t_begin), t_end_(t_end)
   {
     bool bTop_;
-    PbcType pbcType_;
     matrix boxtop_;
     mtop_ = (gmx_mtop_t*)malloc(sizeof(gmx_mtop_t));
     TpxFileHeader header = readTpxHeader(top_path.c_str(), true);
@@ -337,6 +337,7 @@ class CMData
       }
       printf(" - found %s", tmp_mode.c_str());
     }
+    printf("\n");
 
     std::vector<int> num_mol;
     num_mol.push_back(1);
@@ -438,7 +439,7 @@ class CMData
 
     mcut2_ = mol_cutoff_ * mol_cutoff_;
     cut_sig_2_ = (cutoff_ + 0.02) * (cutoff_ + 0.02);
-    if (same_ || cross_) xcm_ = (rvec*)malloc(nindex_ * sizeof(rvec));
+    xcm_ = (rvec*)malloc(nindex_ * sizeof(rvec));
 
     if (same_) frame_same_mat_.resize(natmol2_.size());
     if (intra_ || same_) frame_same_mutex_.resize(natmol2_.size());
@@ -475,7 +476,7 @@ class CMData
     int frnr = 0;
     float progress = 0.0, new_progress = 0.0;
     cmdata::io::print_progress_bar(progress);
-    while (frame_->read_next_frame(trj_) == exdrOK)
+    while (frame_->read_next_frame(trj_, no_pbc_, pbcType_, pbc_) == exdrOK)
     {
       new_progress = static_cast<float>(frnr) / static_cast<float>(frame_->nframe);
       if (new_progress - progress > 0.01)
@@ -503,25 +504,22 @@ class CMData
           #pragma omp parallel for num_threads(std::min(num_threads_, static_cast<int>(frame_cross_mat_[i].size())))
           for ( std::size_t j = 0; j < frame_cross_mat_[i].size(); j++ ) frame_cross_mat_[i][j] = 100.;
         }
-        if (same_ || cross_)
+        #pragma omp parallel for num_threads(std::min(num_threads_, nindex_))
+        for (int i = 0; i < nindex_; i++)
         {
-          #pragma omp parallel for num_threads(std::min(num_threads_, nindex_))
-          for (int i = 0; i < nindex_; i++)
+          clear_rvec(xcm_[i]);
+          double tm = 0.;
+          for (int ii = mols_.block(i).begin(); ii < mols_.block(i).end(); ii++)
           {
-            clear_rvec(xcm_[i]);
-            double tm = 0.;
-            for (int ii = mols_.block(i).begin(); ii < mols_.block(i).end(); ii++)
-            {
-              for (int m = 0; (m < DIM); m++)
-              {
-                xcm_[i][m] += frame_->x[ii][m];
-              }
-              tm += 1.0;
-            }
             for (int m = 0; (m < DIM); m++)
             {
-              xcm_[i][m] /= tm;
+              xcm_[i][m] += frame_->x[ii][m];
             }
+            tm += 1.0;
+          }
+          for (int m = 0; (m < DIM); m++)
+          {
+            xcm_[i][m] /= tm;
           }
         }
         /* start loop for each molecule */
@@ -629,7 +627,7 @@ class CMData
                 end_j_cross = 0;
               }
               if (n_loop_operations_cross == 0) break;
-              if (end_mti_cross == natmol2_.size() - 1) break;
+              if (end_mti_cross == natmol2_.size()) break;
             }
 
             // calculate overhangs and add them

tools/cmdata/src/cmdata/xtc_frame.hpp

@@ -3,6 +3,7 @@
 
 #include "gromacs/math/vec.h"
 #include "gromacs/math/vectypes.h"
+#include <gromacs/pbcutil/pbc.h>
 
 #include <xdrfile_xtc.h>
 
@@ -24,9 +25,17 @@ class Frame {
   Frame(int natom) : natom(natom) { x = (rvec*)malloc(natom * sizeof(rvec)); }
   // ~Frame() { free(x); free(offsets); }
 
-  int read_next_frame(XDRFILE *xd)
+  int read_next_frame(XDRFILE *xd, bool nopbc, PbcType pbc_type, t_pbc *pbc)
   {
     int status = read_xtc(xd, natom, &step, &time, box, x, &prec);
+    if (nopbc)
+    {
+      pbc = nullptr;
+    }
+    else
+    {
+      set_pbc(pbc, pbc_type, box);
+    }
     return status;    
   }
 };

tools/equivalent_atoms/aa_sym.txt

@@ -0,0 +1,36 @@
+# aminoacids
+ARG NH1 NH2
+ASP OD1 OD2
+GLU OE1 OE2
+LEU CD1 CD2
+PHE CD1 CD2
+PHE CE1 CE2
+TYR CD1 CD2
+TYR CE1 CE2
+VAL CG1 CG2
+ALA O1 O2 
+ARG O1 O2 
+ASN O1 O2 
+ASP O1 O2 
+CYS O1 O2 
+GLN O1 O2 
+GLU O1 O2 
+GLY O1 O2 
+HIS O1 O2 
+ILE O1 O2 
+LEU O1 O2 
+LYS O1 O2 
+MET O1 O2 
+PHE O1 O2 
+PRO O1 O2 
+SER O1 O2 
+THR O1 O2 
+TRP O1 O2 
+TYR O1 O2 
+VAL O1 O2 
+# lipids
+POPC C13 C14
+POPC C13 C15
+POPC C14 C15
+POPC O13 O14
+# other as needed

tools/make_mat/make_mat.py

@@ -527,7 +527,8 @@ def calculate_intra_probabilities(args):
 
         if molecule_type == "other":
             # read user pairs
-            user_pairs = [(pair.atom1.idx, pair.atom2.idx, pair.type.epsilon * 4.184) for pair in topology_mego.adjusts]
+            molecule_keys = list(topology_mego.molecules.keys())
+            user_pairs = [(pair.atom1.idx, pair.atom2.idx, pair.type.epsilon * 4.184) for pair in topology_mego.molecules[molecule_keys[i]][0].adjusts]
             user_pairs = [
                 (topology_df[topology_df["mego_ai"] == ai].index[0], topology_df[topology_df["mego_ai"] == aj].index[0], c12)
                 for ai, aj, c12 in user_pairs