from . import constants
def __raw_fepsubst(mol, noncore, cor, non):
"""
Assigns raw fep_subst code to the given molecule. We assign only the following values: 1, 2, 97, 98, and 99, and here we
ignore the other forms of values (i.e., 197, 297, 198, 298, 89, 79, etc.).
:type mol: `structure.Structure`
:param mol: The molecule that we will set the "i_fep_subst" atom property
:type noncore: `set` of `int`
:param noncore: All noncore atoms' indices
:type cor: `list` of `int`
:param cor: All core atoms' indices
:type non: `list` of `int`
:param non: Noncore-but-matched atoms' indices
"""
for e in mol.atom:
e.property[constants.FEP_SUBST] = 1
cor = set(cor)
non = set(non)
for i in noncore:
a = mol.atom[i]
cor_neighbor = set([int(e) for e in a.bonded_atoms]) & cor
if (i in non):
a.property[constants.FEP_SUBST] = 97
else:
a.property[constants.FEP_SUBST] = 99 if (
len(cor_neighbor) > 0) else 2
for j in cor_neighbor:
mol.atom[j].property[constants.FEP_SUBST] = 98
def __is_dummy(atom):
"""
Returns true if this atom can become a dummy atom in an alchemical FEP mutation, or false if not.
"""
fep_subst = atom.property[constants.FEP_SUBST]
return (2 == fep_subst or 99 == fep_subst % 100)
def __fix_97atoms(mol):
"""
We loop through all <97> atoms and check if a <97> atom is connected to more than 1 physical-physical hybrid atoms. If so,
we leave the atom's fep_subst code to be 97, otherwise we change it to 2 (so it becomes a physical-dummy hybrid atom).
We repeat this process until there is no change of the fep_subst code.
"""
should_loop = True
while (should_loop):
should_loop = False
for a in mol.atom:
if (a.property[constants.FEP_SUBST] == 97):
nondummy = [b for b in a.bonded_atoms if (not __is_dummy(b))]
if (len(nondummy) == 1):
a.property[constants.FEP_SUBST] = 99 if (
nondummy[0].property[constants.FEP_SUBST] == 98) else 2
should_loop = True
def __fix_99atoms(atom98, code98, atom99_set):
"""
This function reassigns the 'fep_subst' code for <99> atoms that are bonded to the given <98> atom. Here I use the
<99> nomenclature to mean 99, 199, 299, and so on. -YW
This reassignment is for two purposes:
1. Ensure the code consistency between the two molecules.
2. Adjust some 99 codes to be 89 or 79 or 69.
This reassignment has to be done after that we obtain an overall knowledge of which atoms are <98>, <99>, and <97>.
The assigned 'fep_subst' code is based the given <98> code: if it's 98, <99> will be 99 (or 89, 79, etc.); if it's
198, <99> will be 199 (or 189, 179, etc.), and so on.
:type atom98: `Structure.Atom`
:param atom98: The <98> atom which we will reassign the <99> and <97> atoms that are bonded to
:type code98: `int`
:param code98: The <98> fep_subst code. The value will be one of 98, 198, 298, ... and so on.
:type atom99_set: `set` of `Structure.Atom` objects
:param atom99_set: A set that we use to keep track of the <99> atoms that have already been reassigned. For each atom,
we reassign at most once.
"""
code99 = code98 + 1
atom99 = [
e for e in atom98.bonded_atoms
if (e.property[constants.FEP_SUBST] % 100 == 99)
]
atom99.sort(key=lambda x: x.index)
for e in atom99:
if (e not in atom99_set):
e.property[constants.FEP_SUBST] = code99
atom99_set.add(e)
code99 -= 10
[docs]def get_subst_from_mapping(mol0, mol1):
"""
Given two CTs: `ct0` and `ct1`, with `i_fep_mapping` atom properties properly set for all atoms of `ct1`, sets
the `i_fep_subst` atom properties for both CTs.
"""
# Gets core and non-core matchings.
cor0, cor1 = [], [] # Matched core atoms. List of `int's.
non0, non1 = [], [] # Matched noncore atoms. List of `int's.
for a1 in mol1.atom:
i = a1.property[constants.FEP_MAPPING]
if (i == 0):
continue
a0 = mol0.atom[i]
if (a0.element == a1.element and a0.bond_total == a1.bond_total):
cor0.append(int(a0))
cor1.append(int(a1))
else:
non0.append(int(a0))
non1.append(int(a1))
# Gets core and non-core atoms. Note that some non-core atoms (i.e., the <97> atoms) are matched.
core0 = set(cor0)
core1 = set(cor1)
c1_to_c0 = {c1: c0 for (c0, c1) in zip(cor0, cor1)}
noncore0 = set(range(1, mol0.atom_total + 1))
noncore1 = set(range(1, mol1.atom_total + 1))
noncore0 -= core0
noncore1 -= core1
__raw_fepsubst(mol0, noncore0, cor0, non0)
__raw_fepsubst(mol1, noncore1, cor1, non1)
__fix_97atoms(mol0)
__fix_97atoms(mol1)
atom98 = [(
e.property[constants.FEP_SUBST],
int(e),
) for e in mol1.atom if (e.property[constants.FEP_SUBST] == 98)]
code98 = 98
atom99_set = set()
for junk, i in atom98:
j = c1_to_c0[i]
a1 = mol1.atom[i]
a0 = mol0.atom[j]
a1.property[constants.FEP_SUBST] = code98
a0.property[constants.FEP_SUBST] = code98
__fix_99atoms(a0, code98, atom99_set)
__fix_99atoms(a1, code98, atom99_set)
code98 += 100
# DESMOND-3621
atom97 = [
int(e) for e in mol1.atom if (e.property[constants.FEP_SUBST] == 97)
]
code97 = 97
for i in atom97:
a1 = mol1.atom[i]
j = a1.property[constants.FEP_MAPPING]
a0 = mol0.atom[j]
a1.property[constants.FEP_SUBST] = code97
a0.property[constants.FEP_SUBST] = code97
code97 += 100
# Assigns the i_fep_mapping code.
complete_a1_to_a0 = {a1: a0 for (a0, a1) in zip(cor0 + non0, cor1 + non1)}
for a in mol1.atom:
a.property[constants.FEP_MAPPING] = complete_a1_to_a0[int(a)] if (
a.property[constants.FEP_SUBST] % 100 in [
97,
98,
1,
]) else 0