"""
This program edits a given trajectory as follows:
1. Align each frame to a reference structure as specified via the -ref-mae
option. If this option is not specified, the one from the input CMS file
will be used as the reference structure. The atoms to align can be
specified via the -align-asl option. If it is not specified, the program
will select the atoms automatically based on the following precedences
from high to low:
A. For protein or covalent-ligand FEP systems, the backbone of the mutated
residues will be selected for alignment.
B. If protein is present in the system, all backbone atoms will be selected
for alignment.
C. If it is a FEP system for the solvent leg, the core atoms will be
selected for alignment.
D. If it is a FEP system for the solvent leg, but there are no core atoms
it will be treated as if an absolute binding job and atoms marked
with ABSOLUTE_BINDING_LIGAND will be selected for alignment.
E. All solute atoms will be selected for alignment.
Note that the value of -align-asl will be applied to both the reference
structure and each trajectory frame. If you wish to use a different
atom selection for the reference structure only, you can specify that via
the -ref-asl option.
In any cases, the atoms selected in both the reference structure and each
trajectory frame should match each other. If selected atoms include protein
backbone atoms, all non-backbone atoms will be ignored.
2. For each frame, remove the solvent molecules beyond the region of interest,
which is herein referred to as the 'dew point' and can be specified via
the -dew-asl option. The number of solvent molecules to retain can be
specified via the -n option.
3. In the case of FEP system, atoms corresponding to the other end point will
also be removed.
The default number of solvent molecules retained by this script is 200. This is
approximately the number of water molecules in a spherical droplet of radius
11.5 Angstroms at the temperature of 300 K. To estimate the number of water
molecules to retain, use this table::
nWat: 50 100 200 300 400 500 600 700 800 900 1000
r(A): 7.2 9.1 11.5 13.0 14.3 15.4 16.4 17.3 18.0 18.8 19.6
"""
from typing import List
import numpy
from schrodinger import structure
from schrodinger.application.desmond import cms
from schrodinger.application.desmond import constants
from schrodinger.application.desmond.packages import analysis
from schrodinger.application.desmond.packages import cui
from schrodinger.application.desmond.packages import pfx
from schrodinger.application.desmond.packages import topo
from schrodinger.application.desmond.packages import traj
from schrodinger.application.desmond.packages.timer import Timer
from schrodinger.application.desmond.struc import delete_atom_properties
from schrodinger.structutils import analyze
timer = Timer(threshold=0.01, fmt=" Done. (%.2f sec)", logger=cui.info)
def _is_protein_fep(cms_model):
"""
Return true if the CMS model is for protein FEP simulation, or false if not.
"""
if cms_model.fep_ct is None:
return False
protein_atoms_A = analyze.evaluate_asl(cms_model.ref_ct, "protein")
protein_atoms_B = analyze.evaluate_asl(cms_model.mut_ct, "protein")
return bool(protein_atoms_A) and bool(protein_atoms_B)
def _get_aids_for_alignment(ct, asl) -> List[int]:
"""
This function will be called twice: Once on the reference structure,
the other time on the to-be-aligned structure. Each time, this function
will return a list of AIDs selected by `asl`. We should ensure that the two
results will match each other. This means: (1) the equivalent sets of atoms
are selected, and (2) the atoms are ordered in the same fashion. This
function will guarantee (2), the caller code should check on (1).
Notes:
1. If protein atoms are selected by `asl`, only the backbone atoms' AIDs
are returned, and the AIDs are ordered by connectivity in the peptide
chain (N->C terminus). This is necessary for Protein-FEP, where the
mutated residues in the mutant structure are usually placed at the end
of the atom array.
2. If the `asl` selects both protein and nonprotein atoms, the latter will
be ignored. If you wish to keep the nonprotein atoms, you can call this
function again to only select the nonprotein atoms, and then combine
the result with the previous one.
"""
bb_aids = analyze.evaluate_asl(
ct, f'({asl}) and backbone and protein and not a.e H')
if bb_aids:
PROP_ORIG_INDEX = 'i_des_tmp_original_index'
for a in ct.atom:
a.property[PROP_ORIG_INDEX] = a.index
bb_ct = ct.extract(bb_aids)
ordered_bb_atoms = []
for res in structure.get_residues_by_connectivity(bb_ct):
# `.pdbres` is a 4-char string. If it's a standard PDB residue
# name (which strictly has 3 char's), then `.pdbres`'s value is
# the standard name plus a trailing space.
if res.pdbres not in ['ACE ', 'NMA ']:
ordered_bb_atoms.extend([
res.getBackboneNitrogen(),
res.getAlphaCarbon(),
res.getCarbonylCarbon(),
res.getBackboneOxygen()
])
result = [a.property[PROP_ORIG_INDEX] for a in ordered_bb_atoms if a]
delete_atom_properties(ct, PROP_ORIG_INDEX)
return result
return analyze.evaluate_asl(ct, asl)
def _fix_res_number(ct):
# Resets the residue numbers to ensure residue numbers start with 1
# and are contiguous. The only reason that they are not so now is
# because they were like that in the input <cms> file.
resnum = 1
for mol in ct.molecule:
# Atoms belonging to the same molecule may NOT be contiguous in
# the atom array. We assume each solvent molecule has only 1
# "residue", which is usually true.
for atom in mol.atom:
atom.resnum = resnum
resnum += 1
def _delete_extra_solvent(ct, num_solvent):
# Finds out the number of pseudoatoms of each solvent molecule.
solvent_atom_total = ct.atom_total
solvent_pseudo_total = len(ct.ffio.pseudo)
solvent_molecule_total = len(ct.molecule)
num_atoms_per_solvent = solvent_atom_total // solvent_molecule_total
num_pseudo_per_solvent = solvent_pseudo_total // solvent_molecule_total
assert num_pseudo_per_solvent * solvent_molecule_total == \
solvent_pseudo_total, \
f"""num_pseudo_per_solvent * solvent_molecule_total == solvent_pseudo_total:
number of pseudoatoms per solvent molecule: {num_pseudo_per_solvent}
total number of solvent molecules: {solvent_molecule_total}
total number of pseudoatoms: {solvent_pseudo_total}"""
ct.deleteAtoms(
range(num_atoms_per_solvent * num_solvent + 1, solvent_atom_total + 1))
for i in range(solvent_pseudo_total, num_pseudo_per_solvent * num_solvent,
-1):
del ct.ffio.pseudo[i]
_fix_res_number(ct)
def _reduce_solvent(msys_model, cms_model, fep_lambda, num_solvent, tr):
"""
We divide the overall procedure into 3 steps:
1. Recenter the trajectory to find out solvent molecules close to the
origin.
2. Align the trajectory.
3. Reduce the solvent by keeping the closest solvent molecules already
selected at Step 1.
Note that recentering at Step 1 and alignment at Step 2 both can alter the
3D coordinates in the trajectory such that they will mess up each other's
result. Although Step 2 has nothing to do with reducing the solvent, we
have to let it happen between Steps 1 and 3. This is better implemented as
a coroutine where we suspend this function after Step 1 and resume it after
Step 2 that will be done at the caller code level.
This coroutine is in the form of Python generator. Use example::
reduce_solvent_by_2_steps = _reduce_solvent(...)
next(reduce_solvent_by_2_steps) # Step 1
# Step 2: Align the trajectory.
next(reduce_solvent_by_2_steps) # Step 3
The atoms that we will trim off are solvent molecules that are far away
from the "dew point", plus (in the case of FEP) the alchemical molecule in
the other lambda state (which is 0 if `fep_lambda` = 1, or 1 if
`fep_lambda` = 0).
Handlings of corner cases:
- If there are less solvent molecules than the requested number, this
function will simply keep all (i.e., as many as possible) solvent
molecules.
- If there is no solvent CT at all, this function will return the input
`tr` and print out a warning message saying there is no solvent CT.
:type fep_lambda: 0, 1, or None
:param fep_lambda: Lambda window (either 0 or 1) if `cms_model` is a FEP
system. The value is `None` if `cms_model` is a non-FEP system.
:param num_solvent: Number of solvent molecules to keep.
:rtype: (cms.Cms, List(traj.Frame))
:return: The CMS model and the trajectory of the parched system
"""
# Notes:
# 1. We only consider molecules from the solvent CT(s) as the solvent
# molecules. There might be the same types of molecules in other
# component CTs, we won't count them as solvent molecules.
# 2. For mixed solvent, we may have more than 1 solvent CTs. But this
# function is not full-fledged to handle mixed solvent.
solvent_cts = [
ct for ct in cms_model.comp_ct
if ct.property[constants.CT_TYPE] == constants.CT_TYPE.VAL.SOLVENT
]
# FIXME: For simplicity, from this point on, we assume there is only
# one solvent CT. We need extra fix for mixed solvents.
# -YW 5/23/2018
if not solvent_cts:
cui.warn("Solvent CT not found. No solvent molecules to parch.")
else:
# Ensures to not select inactive atoms.
solvent_aids = cms_model.select_atom("solvent")
first_solvent_aid = min(solvent_aids)
first_solvent_molecule = cms_model.atom[first_solvent_aid].getMolecule()
solvent_molecule_total = sum(len(e.molecule) for e in solvent_cts)
if num_solvent < 1:
cui.info(" Will strip off all solvent molecules...")
else:
cui.info(" Calculating distances of all solvent molecules "
"to the dew point...")
heavy_atoms_per_solvent = [
a for a in first_solvent_molecule.atom if a.atomic_number > 1
]
if len(heavy_atoms_per_solvent) == 1:
# For the case where each solvent molecule has only 1 heavy
# atom, we can use `AtomicPosition` analyzer, which is ~2x
# faster than `MoleculeWiseCom`.
ana = analysis.AtomicPosition(msys_model,
cms_model,
asl="solvent_heavy_atom")
else:
ana = analysis.MoleculeWiseCom(msys_model, cms_model, "solvent")
def progress_report(i, *_):
print(f" processing frame# {i}...")
ana_result = analysis.analyze(tr,
ana,
progress_feedback=progress_report)
# Because the simulation system is already centered at the "dew"
# point, we just need to calculate the r^2, namely the squared
# Euclidean-norm of the position vectors, of the COMs.
# FIXME: This is a bit corner-cutting, ignoring the PBC. But this
# has been used for a while. So far there is no sign that
# this is a real problem in practice. If it turns out to be
# so, we should fix. -YW 5/22/2018
tr_r2 = [(sol_pos**2).sum(axis=-1) for sol_pos in ana_result]
# Do NOT reduce the solvent yet.
# Temporarily suspends this function until we have aligned all the frames.
yield
# Creates a new CMS model.
new_cms_model = cms_model.copy()
if cms_model.need_msys:
_ORIG_AID_PROP = 'i_des_orig_aid'
# Records the full system CT's AIDs on each component CT's atoms. This is
# to map the subsystem AIDs to original-system AIDs.
for a0, a1 in zip(topo.comp_atoms(new_cms_model), new_cms_model.atom):
a0.property[_ORIG_AID_PROP] = a1.index
alchem_ct_to_skip = None
if fep_lambda is not None:
# Finds out the alchemical CT to *delete*.
all_aids = cms_model.select_atom("all")
alchem_ct = (cms_model.ref_ct if fep_lambda == 1 else cms_model.mut_ct)
alchem_ct_to_skip = cms_model.comp_ct.index(alchem_ct)
alchem_aids = cms_model.get_fullsys_ct_atom_index_range(
alchem_ct_to_skip)
if cms_model.need_msys:
new_cms_model, _ = topo.extract_subsystem(
new_cms_model,
"(a.n %s)" % " ".join(map(str,
set(all_aids) - set(alchem_aids))))
# Finds the solvent CT.
solvent_ct_idx = None
if solvent_cts:
# FIXME: We assume there is only one solvent CT.
for index_solvent_ct, ct in enumerate(new_cms_model.comp_ct):
if ct.property[constants.CT_TYPE] == constants.CT_TYPE.VAL.SOLVENT:
solvent_ct_idx = index_solvent_ct
break
if cms_model.need_msys:
if num_solvent < 1:
del new_cms_model.comp_ct[index_solvent_ct]
else:
_delete_extra_solvent(ct, num_solvent)
orig_active_total = new_cms_model.active_total
new_cms_model.active_total = new_cms_model.comp_atom_total
# This is to avoid syncing twice since `active_total` setter may sync.
# It also triggers sync when comp_atom_total after parch happens to be
# equal to active_total (it is the reverse of the Cms.active_total
# setter logic).
if (new_cms_model._show_inactive or
orig_active_total == new_cms_model.comp_atom_total):
new_cms_model.synchronize_fsys_ct()
def _extract_cms_with_some_solvent(cms_model, solvent_ct_idx,
num_solvent, alchem_ct_to_skip):
new_solvent_ct_idx = None
comp_sel = []
tmp_id = 0
for ct_idx, ct in enumerate(cms_model.comp_ct):
if ct_idx == alchem_ct_to_skip:
comp_sel.append([])
continue
if ct_idx == solvent_ct_idx:
solvent_selection = list(
range(1,
len(ct.molecule[1]) * num_solvent + 1))
if len(solvent_selection):
new_solvent_ct_idx = tmp_id
comp_sel.append(solvent_selection)
else:
comp_sel.append(list(range(1, ct.atom_total + 1)))
tmp_id += 1
new_cms_model, kept_gids = cms_model.extract_with_comp_selection(
comp_sel)
if new_solvent_ct_idx is not None:
_fix_res_number(new_cms_model.comp_ct[new_solvent_ct_idx])
new_cms_model.synchronize_fsys_ct()
return new_cms_model, kept_gids, new_solvent_ct_idx
if cms_model.need_msys:
new_msys_model, new_cms_model = topo.read_cms(
from_string=new_cms_model.write_to_string())
else:
new_cms_model, kept_gids, new_solvent_ct_idx = _extract_cms_with_some_solvent(
cms_model, solvent_ct_idx, num_solvent, alchem_ct_to_skip)
# Creates the new trajectory.
if cms_model.need_msys:
kept_gids = topo.matched_gids(cms_model, new_cms_model, new_msys_model)
if solvent_cts:
if cms_model.need_msys:
cui.info(" Deleting solvent molecules for each frame...")
# Removes solvent GIDs from `kept_gids`. We will figure out the solvent
# GIDs for each frame.
solvent_gids = topo.asl2gids(cms_model, "solvent")
first_solvent_gid = min(solvent_gids)
last_solvent_gid = max(solvent_gids)
kept_gids = [
i for i in kept_gids
if i < first_solvent_gid or i > last_solvent_gid
]
new_tr = []
for i, fr in enumerate(tr):
print(f" processing frame# {i}...")
if (not solvent_cts) or (num_solvent < 1):
solvent_kept_gids = []
new_tr.append(fr.reduce(kept_gids))
else:
# Partitions `fr_r2` on value, and then gets the molecule indices of the solvent
# molecules that have the smallest `num_solvent` r^2 values.
fr_r2 = tr_r2[i]
kept_solvent_molecule_indices = numpy.argpartition(
fr_r2, num_solvent)[:num_solvent]
cms_model.set_nactive_gids(fr.nactive, fr.natoms)
solvent_asl = "m.n %s" % " ".join(
map(
str, kept_solvent_molecule_indices +
first_solvent_molecule.number))
if cms_model.need_msys:
solvent_kept_gids = topo.asl2gids(cms_model,
solvent_asl,
include_pseudoatoms=True)
new_tr.append(fr.reduce(kept_gids + sorted(solvent_kept_gids)))
else:
offset = new_cms_model.id_maps[new_solvent_ct_idx].start_gid
aids = cms_model.select_atom(solvent_asl)
for sol_index, orig_gid in enumerate(
cms_model.convert_to_gids(aids, with_pseudo=True),
offset):
kept_gids[sol_index] = orig_gid
new_tr.append(fr.reduce(kept_gids))
if cms_model.need_msys:
delete_atom_properties([new_cms_model] + new_cms_model.comp_ct,
_ORIG_AID_PROP)
yield new_cms_model, new_tr
_ATOM_PROPS_TO_REMOVE = {
'r_ffio_x_vel',
'r_ffio_y_vel',
'r_ffio_z_vel',
constants.ORIGINAL_INDEX,
constants.REST_PROPERTIES.COMPLEX_HOTREGION,
constants.REST_PROPERTIES.SOLVENT_HOTREGION,
constants.FEP_MAPPING,
}.union(set(constants.REFERENCE_COORD_PROPERTIES))
_CT_PROPS_TO_REMOVE = [
constants.MOLTYPE, constants.FEPIO_STAGE, constants.FEP_FRAGNAME
]
# FIXME: Should these functions (clear and reset) be placed in mmshare?
def _clear_fep_props_for_ct(ct: structure.Structure) -> None:
for p in _CT_PROPS_TO_REMOVE:
if p in ct.property:
del ct.property[p]
for atom in ct.atom:
for p in _ATOM_PROPS_TO_REMOVE:
if p in atom.property:
del atom.property[p]
def _clear_properties(cms_model: cms.Cms) -> None:
for ct in [cms_model.fsys_ct, cms_model._raw_fsys_ct] + cms_model.comp_ct:
_clear_fep_props_for_ct(ct)
def _get_cml():
cml = cui.CommandLine(
(cui.REQUIRE_CMS_ONLY + cui.REQUIRE_TRJ + cui.REQUIRE_OUT +
cui.REF_MAE + cui.SLICE_TRJ + cui.OUT_TRJ_FORMAT),
description=__doc__)
cml.edit_argument_help(
flag="-ref-mae",
new_help=" If not specified, the default reference structure will be "
"that from the <cms> argument.",
edit="append")
cml.add_argument(
"-ref-asl",
metavar='<ASL>',
default=None,
type=cui.validate_and_enclose_asl,
help="ASL expression to select atoms in the reference structure. "
"If not specified, the same ASL expression as that of -align-asl "
"will be evaluated against the reference structure to get the "
"reference coordinates.")
cml.add_argument(
"-align-asl",
metavar='<ASL>',
default=None,
type=cui.validate_and_enclose_asl,
help="ASL expression to select atoms in each trajectory frame. The "
"atoms selected in each frame will be aligned to atoms in the "
"reference structure selected by -ref-asl. If not specified, the "
"program will select atoms automatically (see above for detail).")
cml.add_argument(
"-dew-asl",
metavar='<ASL>',
required=True,
type=cui.validate_and_enclose_asl,
help="ASL expression to specify the atoms around which a solvent "
"droplet ('dew') will be retained. (Required)")
cml.add_argument('-nwaters',
metavar='<int>',
type=int,
help="(DEPRECATED! Please use -n instead.)")
cml.add_argument(
'-n',
metavar='<int>',
dest='num_solvent',
type=int,
default=None, # FIXME: s/None/200, when -nwaters is deleted.
help="Number of solvent molecules to keep. Default: 200")
cml.add_argument("-fep-lambda",
type=int,
default=None,
choices=[0, 1],
help="Specify the FEP lambda end point to keep.")
cml.add_argument(
'-protein-fep',
action='store_true',
default=False,
help="(DEPRECATED! This option has no effects now. "
"The program will automatically detect if it is a protein FEP "
"system or not.)")
cml.edit_argument_help(
flag="-output-trajectory-format",
new_help=(", which means the same format as the input trajectory."),
edit="append")
return cml
def _main_impl(argv=None):
cml = _get_cml()
arg = cml.parse_args(argv)
# Sanity checks on command line arguments
# Ensures some options' defaults are set correctly.
cms_model, cms_fname = arg.cms
msys_model = None
if cms_model.need_msys:
msys_model, cms_model = topo.read_cms(cms_fname)
tr, tr_fname = arg.trj
if arg.out_trj_format == traj.Fmt.AUTO:
arg.out_trj_format = traj.get_fmt(tr_fname)
if arg.fep_lambda is None and cms_model.fep_ct:
cui.error("The input system contains an FEP block, please use "
"-fep-lambda to specify an end-point of interest.")
if arg.ref_mae is None:
arg.ref_mae = cms_model.fsys_ct, cms_fname
if arg.nwaters is not None:
cui.warn("-nwaters has been DEPRECATED. Please use -n instead.")
if arg.num_solvent is None:
arg.num_solvent = arg.nwaters or 200
if arg.align_asl is None:
if arg.fep_lambda is not None:
# FEP system
if _is_protein_fep(cms_model):
# For protein FEP systems, alignment is on mutated residue(s).
ct = (cms_model.ref_ct
if arg.fep_lambda == 0 else cms_model.mut_ct)
ct_index = cms_model.comp_ct.index(ct)
aid_range = cms_model.get_fullsys_ct_atom_index_range(ct_index)
arg.align_asl = "(a.n %d-%d)" % (aid_range[0], aid_range[-1])
else:
# Small compound FEP systems
if cms_model.select_atom("protein"):
# If it's a complex system, alignment is on the protein.
arg.align_asl = "protein"
else:
# If it's a solvent system, alignment is on the core atoms.
# If there are no core atoms, this is an absolute binding
# calculation, so use the ligand instead
arg.align_asl = f"atom.{constants.FEP_MAPPING} > 0 or a.{constants.FEP_ABSOLUTE_BINDING_LIGAND} 1"
else:
arg.align_asl = f"atom.{constants.CT_TYPE} '{constants.CT_TYPE.VAL.SOLUTE}'"
cui.info(f'Automatic value for -align-asl: "{arg.align_asl}"')
if arg.ref_asl is None:
arg.ref_asl = arg.align_asl
cui.info(f"Automatic value for -ref-asl: {arg.ref_asl}")
# By now, all arguments should have correct values.
with timer(f"1. Centering the entire system on atoms: {arg.dew_asl}..."):
dew_aids = cms_model.select_atom(arg.dew_asl)
cui.info(" %d atoms selected by ASL expression: %s" %
(len(dew_aids), arg.dew_asl))
if not dew_aids:
cui.error(f"No atoms selected by -dew-asl: {arg.dew_asl}.")
if cms_model.need_msys:
tr = topo.center(msys_model, topo.aids2gids(cms_model, dew_aids),
tr)
else:
tr = pfx.center(cms_model.glued_topology,
cms_model.convert_to_gids(dew_aids), tr)
reduce_solvent_by_2_steps = _reduce_solvent(msys_model, cms_model,
arg.fep_lambda, arg.num_solvent,
tr)
with timer(f"2. Selecting closest {arg.num_solvent} solvent "
f"molecules to atoms: {arg.dew_asl}..."):
next(reduce_solvent_by_2_steps)
with timer("3. Selecting atoms in reference structure for alignment..."):
ref_mae, _ = arg.ref_mae
ref_align_aids = _get_aids_for_alignment(ref_mae, arg.ref_asl)
num_ref_align_aids = len(ref_align_aids)
if 3 > num_ref_align_aids:
cui.error(f"{num_ref_align_aids} atoms selected in reference "
f"structure by '-ref-asl \"{arg.ref_asl}\"'. At least 3 "
f"atoms should be selected.")
cui.info(" %d atoms were selected." % num_ref_align_aids)
with timer("4. Selecting atoms in <cms> structure for alignment..."):
align_aids = _get_aids_for_alignment(cms_model.fsys_ct, arg.align_asl)
num_align_aids = len(align_aids)
if 3 > num_align_aids:
cui.error(f"{num_align_aids} atoms selected in <cms> structure "
f"by '-align-asl \"{arg.align_asl}\"'. At least 3 atoms "
f"should be selected.")
cui.info(f" {num_align_aids} atoms were selected.")
# Checks if `ref_align_aids` matches `align_aids`.
if num_ref_align_aids != num_align_aids:
cui.error("The number of atoms selected by the -align-asl in reference "
"structure does not match that in <cms> structure.")
for i, j in zip(align_aids, ref_align_aids):
a = cms_model.atom[i].atomic_weight
b = ref_mae.atom[j].atomic_weight
if a != b:
cui.warn(f"Aligned atoms are of different types: "
f"Atomic weights: {a} vs {b}")
with timer("5. Aligning trajectory on reference structure..."):
# NOTE: getXYZ() return a numpy array with 0-based index
ref_align_pos = ref_mae.getXYZ()[numpy.array(ref_align_aids) - 1]
align_weights = [ref_mae.atom[i].atomic_weight for i in ref_align_aids]
if cms_model.need_msys:
align_gids = topo.aids2gids(cms_model,
align_aids,
include_pseudoatoms=False)
tr = topo.superimpose(msys_model,
align_gids,
tr,
ref_pos=ref_align_pos,
weights=align_weights)
else:
align_gids = cms_model.convert_to_gids(align_aids,
with_pseudo=False)
tr = pfx.superimpose(cms_model.glued_topology,
align_gids,
tr,
ref_pos=ref_align_pos,
weights=align_weights)
with timer("6. Extracting solute and the solvent "
"molecules selected at step 2..."):
out_cms_model, out_tr = next(reduce_solvent_by_2_steps)
if cms_model.need_msys:
err_msg = topo.check_consistency(out_cms_model, out_tr[-1])
if err_msg:
cui.error(
f"Found inconsistency between the parched CMS and the "
f"parched trajectory:\n{err_msg}.")
topo.update_cms(out_cms_model, out_tr[-1], update_pseudoatoms=True)
else:
out_cms_model.update_with_frame(out_tr[-1])
cui.info(" %d atoms now in parched system (%d atoms originally)" %
(out_cms_model.atom_total, cms_model.atom_total))
with timer("7. Writing the new CMS file and trajectory..."):
tr_extname = traj.get_trajectory_extname(arg.out_trj_format, tr_fname)
out_trj_fname = f"{arg.out}{tr_extname}"
out_cms_fname = f"{arg.out}-out.cms"
traj.write_traj(out_tr, out_trj_fname, format=arg.out_trj_format)
_clear_properties(out_cms_model)
out_cms_model.fix_filenames(out_cms_fname, out_trj_fname)
out_cms_model.write(out_cms_fname)
cui.info(f" Output files:\n {out_cms_fname}\n {out_trj_fname}")
return out_cms_fname, out_trj_fname
[docs]def main(argv=None):
with timer(fmt="All done. (%.2f sec)"):
return _main_impl(argv)