fix ENH issue #454 qm atom mask selection

Valdes-Tresanco-MS · Jan 17, 2024 · 027a19b · 027a19b
1 parent 71a07d1
commit 027a19b
Show file tree

Hide file tree

Showing 4 changed files with 97 additions and 57 deletions.
diff --git a/GMXMMPBSA/createinput.py b/GMXMMPBSA/createinput.py
@@ -199,6 +199,12 @@ def create_inputs(INPUT, prmtop_system, pre):
                 lig_input = deepcopy(INPUT)
                 (com_input['gb']['qmmask'], rec_input['gb']['qmmask'],
                  lig_input['gb']['qmmask']) = prmtop_system.Mask(INPUT['gb']['qm_residues'], in_complex=False)
+
+                if INPUT['gb']['exclude_backbone']:
+                    com_input['gb']['qmmask'] = f"({com_input['gb']['qmmask']} & !@N,CA,C,O)"
+                    rec_input['gb']['qmmask'] = f"({rec_input['gb']['qmmask']} & !@N,CA,C,O)"
+                    lig_input['gb']['qmmask'] = f"({lig_input['gb']['qmmask']} & !@N,CA,C,O)"
+
                 if not com_input['gb']['qmmask']:
                     raise AmberError('No valid QM residues chosen!')
                 com_input['gb']['qm_theory'] = "'%s'" % com_input['gb']['qm_theory']

diff --git a/GMXMMPBSA/input_parser.py b/GMXMMPBSA/input_parser.py
@@ -450,6 +450,8 @@ def Parse(self, filename):
                            ['qm_theory', str, '', 'Semi-empirical QM theory to use'],
                            ['qm_residues', str, '', 'Residues to treat with QM'],
 
+                           ['exclude_backbone', int, 0, 'Exclude backbone atoms from residues to treat with QM'],
+
                            # TODO: deprecated since 1.5.0. Automatic charge assignment
                            ['qmcharge_com', int, 0, 'Charge of QM region in complex'],
                            ['qmcharge_lig', int, 0, 'Charge of QM region in ligand'],

diff --git a/GMXMMPBSA/main.py b/GMXMMPBSA/main.py
@@ -1017,6 +1017,8 @@ def check_for_bad_input(self, INPUT=None):
         if INPUT['gb']['ifqnt'] == 0 and (INPUT['gb']['qm_theory'] or INPUT['gb']['qm_residues']):
             logging.warning('qm_theory/qm_residues variable has been defined, however the potential function is '
                             'strictly classical (ifqnt=0). Please, set ifqnt=1 if you want to use Use QM/MM')
+        if INPUT['gb']['exclude_backbone'] not in [0, 1]:
+            GMXMMPBSA_ERROR('exclude_backbone must be 0 or 1!', InputError)
         if (
                 not INPUT['gb']['molsurf']
                 and (INPUT['gb']['msoffset'] != 0 or INPUT['gb']['probe'] != 1.4)