#!/usr/bin/env python
#
# Copyright (c) 2024 Authors and contributors
# (see the AUTHORS.rst file for the full list of names)
#
# Released under the GNU Public Licence, v3 or any higher version
# SPDX-License-Identifier: GPL-3.0-or-later
"""Build universes from template molecules."""
import logging
from typing import Optional
import MDAnalysis as mda
import numpy as np
from tqdm import tqdm
from .models import empty
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)
def _renumber_projectile_resids(
SolvatedUniverse: mda.Universe, nAtomsTarget: int
) -> mda.Universe:
"""Renumber residues after the target so resids are contiguous and monotonic.
Target residues (the first `nAtomsTarget` atoms) keep their original resids.
Residues from the projectile atoms are renumbered starting at
``target.resids[-1] + 1`` (or ``1`` when the target is empty), so the
returned universe has no gaps or duplicates among the projectile resids.
"""
n_total_res = len(SolvatedUniverse.residues)
if nAtomsTarget == 0:
start = 1
n_target_res = 0
else:
target = SolvatedUniverse.atoms[:nAtomsTarget]
# Use max() rather than [-1] so we don't collide with an
# out-of-order target resid (e.g. user-supplied [5, 2, 3]).
start = target.residues.resids.max() + 1
n_target_res = len(target.residues)
n_proj_res = n_total_res - n_target_res
if n_proj_res > 0:
projectile = SolvatedUniverse.atoms[nAtomsTarget:]
projectile.residues.resids = np.arange(start, start + n_proj_res)
return SolvatedUniverse
def tile_universe(
universe: mda.Universe,
n_x: int,
n_y: int,
n_z: int,
) -> mda.Universe:
"""Returns a new Universe with `n_x * n_y * n_z` copies of the input."""
box = universe.dimensions[:3]
copied = []
i = 0
for x in tqdm(range(n_x)):
for y in range(n_y):
for z in range(n_z):
u_ = universe.copy()
move_by = box * (x, y, z)
u_.residues.resids += len(universe.residues) * i
u_.atoms.translate(move_by)
copied.append(u_.atoms)
i += 1
new_universe = mda.Merge(*copied)
new_box = box * (n_x, n_y, n_z)
new_universe.dimensions = list(new_box) + [90] * 3
return new_universe
def pos_random(InsertionDomain: np.ndarray) -> np.ndarray:
"""Returns a random position within the given domain."""
return np.array(
np.random.rand(3) * (InsertionDomain[3:6] - InsertionDomain[0:3])
+ InsertionDomain[0:3],
dtype=np.float32,
)
def rot_random() -> tuple:
"""Returns a random rotation angle and vector, sampled uniformly on a sphere."""
u_1, u_2, u_3 = np.random.rand(3)
theta, phi = np.arccos(2 * u_1 - 1), 2 * np.pi * u_2
rot_vec = np.array(
[np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)]
)
rot_angle = 360 * u_3
return rot_angle, rot_vec
[docs]
def SolvateCylinder(
TargetUniverse: mda.Universe,
ProjectileUniverse: mda.Universe,
n: int = 1,
density: Optional[float] = None,
pos: Optional[np.ndarray] = None,
radius: Optional[float] = None,
min: float = 0,
max: Optional[float] = None,
dim: int = 2,
distance: float = 1.25,
tries: int = 1000,
fudge_factor: float = 1,
) -> mda.Universe:
"""Fill a cylindrical region of a target with copies of a projectile.
Internally builds a small saturated patch of projectiles, tiles it across
the cylinder's bounding box, then prunes any copies that fall outside the
cylinder or overlap atoms in ``TargetUniverse``. For most use cases this is
orders of magnitude faster than :func:`InsertCylinder`.
Parameters
----------
TargetUniverse : MDAnalysis.core.universe.Universe
Universe to solvate. May be empty; in that case
``TargetUniverse.dimensions`` still defines the simulation cell of the
returned universe.
ProjectileUniverse : MDAnalysis.core.universe.Universe
Molecule that is inserted repeatedly.
n : int, default 1
Number of projectile copies to insert. Ignored when ``density`` is
given.
density : float, optional
Target number density of projectiles, in molecules / ų. When set,
``n`` is computed from the cylinder volume and ``n`` is ignored.
pos : array_like of shape (3,), optional
Centre of the cylinder, in Å. Defaults to the centre of geometry of
``TargetUniverse`` (or the centre of its box if the target is empty).
The component along ``dim`` is overridden by ``min``.
radius : float, optional
Cylinder radius, in Å. Defaults to half of the smallest box edge of
``TargetUniverse``.
min, max : float, optional
Lower and upper bound of the cylinder along axis ``dim``, in Å.
``max`` defaults to ``TargetUniverse.dimensions[dim]``.
dim : {0, 1, 2}, default 2
Index of the axis along which the cylinder extends (0 = x, 1 = y,
2 = z).
distance : float, default 1.25
Minimum allowed distance (Å) between an inserted projectile and any
atom of the target.
tries : int, default 1000
Maximum number of random placement attempts used when building the
seed patch.
fudge_factor : float, default 1.0
Multiplier on the number of projectiles packed into the seed patch.
Increased automatically and the function recurses when too few
projectiles survive overlap pruning.
Returns
-------
MDAnalysis.core.universe.Universe
New universe containing the original target atoms followed by the
inserted projectile copies.
See Also
--------
InsertCylinder : Slower per-molecule variant with full placement control.
SolvatePlanar : Equivalent for rectangular regions.
"""
logger.info(f"The fudge factor is {fudge_factor}")
if max is None:
max = TargetUniverse.dimensions[dim]
nAtomsProjectile = ProjectileUniverse.atoms.n_atoms
dimensionsTarget = TargetUniverse.dimensions.copy()
if pos is None:
if TargetUniverse.atoms.n_atoms == 0:
pos = dimensionsTarget[:3] / 2
else:
pos = TargetUniverse.atoms.center_of_geometry()
pos[dim] = min
if radius is None:
radius = np.min(dimensionsTarget) / 2
if density is not None:
n = np.floor(density * (2 * radius) ** 2 * (max - min))
solvate_by_density_flag = True
else:
solvate_by_density_flag = False
dimensionsTarget = TargetUniverse.dimensions.copy()
nAtomsTarget = TargetUniverse.atoms.n_atoms
nAtomsProjectile = ProjectileUniverse.atoms.n_atoms
InsertionDomain = np.array(
[pos[0] - radius, pos[1] - radius, min, pos[0] + radius, pos[1] + radius, max],
dtype=np.float32,
)
InsertionVolume = (max - min) * np.pi * radius**2
density = n / InsertionVolume
SolvatedUniverse = SolvatePlanar(
TargetUniverse,
ProjectileUniverse,
0,
density,
xmin=InsertionDomain[0],
ymin=InsertionDomain[1],
zmin=InsertionDomain[2],
xmax=InsertionDomain[3],
ymax=InsertionDomain[4],
zmax=InsertionDomain[5],
distance=distance,
tries=tries,
fudge_factor=fudge_factor,
)
dims = SolvatedUniverse.dimensions
TargetAtoms = SolvatedUniverse.atoms[:nAtomsTarget]
ProjectileAtoms = SolvatedUniverse.atoms[nAtomsTarget:]
atomsInside = (
np.linalg.norm((ProjectileAtoms.positions - pos)[:, :2], axis=1) < radius
)
if TargetAtoms.n_atoms == 0:
SolvatedUniverse = ProjectileAtoms[atomsInside].residues.atoms
else:
SolvatedUniverse = mda.Merge(
TargetAtoms, ProjectileAtoms[atomsInside].residues.atoms
)
SolvatedUniverse.dimensions = dims
logger.info("Resulting number of atoms:", SolvatedUniverse.atoms.n_atoms)
logger.info(
"Resulting number of projectiles:",
(SolvatedUniverse.atoms.n_atoms - nAtomsTarget) / nAtomsProjectile,
)
missingProjectiles = int(
((n * nAtomsProjectile + nAtomsTarget) - SolvatedUniverse.atoms.n_atoms)
/ nAtomsProjectile
)
logger.info("Missing", missingProjectiles, "Projectiles.")
if solvate_by_density_flag:
logger.info(
f" {SolvatedUniverse.atoms.n_atoms - nAtomsTarget} projectiles inserted"
)
return _renumber_projectile_resids(SolvatedUniverse, nAtomsTarget)
if missingProjectiles > 0:
logger.info("Missing", missingProjectiles, "Projectiles.")
logger.info("Adjusting fudge factor and trying again.")
new_fudge_factor = fudge_factor + 0.5
return SolvateCylinder(
TargetUniverse,
ProjectileUniverse,
n,
density=None,
pos=pos,
radius=radius,
min=min,
max=max,
dim=dim,
distance=distance,
tries=tries,
fudge_factor=new_fudge_factor,
)
if missingProjectiles < 0:
nonTargetAtoms = SolvatedUniverse.atoms[nAtomsTarget:]
logger.info("Too many projectiles inserted:", -missingProjectiles)
logger.info(nonTargetAtoms.n_atoms)
logger.info(nonTargetAtoms.residues.n_residues)
logger.info(np.unique(nonTargetAtoms.residues.resids).shape)
logger.info("Removing", -missingProjectiles, "randomly selected projectiles.")
ToBeRemoved = nonTargetAtoms.residues[
np.random.choice(
np.arange(len(nonTargetAtoms.residues)),
-missingProjectiles,
replace=False,
)
]
SolvatedUniverse = mda.Merge(SolvatedUniverse.atoms - ToBeRemoved.atoms)
SolvatedUniverse.dimensions = dimensionsTarget
logger.info("Final number of atoms:", SolvatedUniverse.atoms.n_atoms)
return _renumber_projectile_resids(SolvatedUniverse, nAtomsTarget)
logger.info("All projectiles inserted correctly")
return _renumber_projectile_resids(SolvatedUniverse, nAtomsTarget)
[docs]
def SolvatePlanar(
TargetUniverse: mda.Universe,
ProjectileUniverse: mda.Universe,
n: int = 1,
density: Optional[float] = None,
xmin: int = 0,
ymin: int = 0,
zmin: int = 0,
xmax: Optional[float] = None,
ymax: Optional[float] = None,
zmax: Optional[float] = None,
distance: float = 1.25,
solvate_factor: int = 100,
fudge_factor: float = 1.0,
tries: int = 1000,
) -> mda.Universe:
"""Fill a rectangular region of a target with copies of a projectile.
Internally builds a small saturated patch of projectiles, tiles it across
the insertion box, then prunes any copies that overlap atoms in
``TargetUniverse``. This is orders of magnitude faster than
:func:`InsertPlanar` for large solvent counts and is the recommended way
to solvate a target with thousands of solvent molecules.
Parameters
----------
TargetUniverse : MDAnalysis.core.universe.Universe
Universe to solvate. May be empty; ``TargetUniverse.dimensions``
defines the simulation cell of the returned universe.
ProjectileUniverse : MDAnalysis.core.universe.Universe
Molecule that is inserted repeatedly.
n : int, default 1
Number of projectile copies to insert. Ignored when ``density`` is
given.
density : float, optional
Target number density of projectiles, in molecules / ų. When set,
``n`` is computed from the volume of the insertion box and ``n`` is
ignored.
xmin, ymin, zmin : float, default 0
Lower bounds of the insertion box, in Å.
xmax, ymax, zmax : float, optional
Upper bounds of the insertion box, in Å. Each defaults to the
corresponding component of ``TargetUniverse.dimensions``.
distance : float, default 1.25
Minimum allowed distance (Å) between an inserted projectile and any
atom of the target. Tile copies closer than ``distance`` are removed
after tiling.
solvate_factor : int, default 100
Target number of projectiles in each tiled sub-box. Larger values
reduce the number of tiles and the cost of the per-tile saturation
step; smaller values reduce peak memory.
fudge_factor : float, default 1.0
Multiplier on ``solvate_factor`` controlling how aggressively the
seed patch is packed. Increased automatically and the function
recurses when too few projectiles survive overlap pruning.
tries : int, default 1000
Base number of random placement attempts used when packing the seed
patch (internally scaled by 1000).
Returns
-------
MDAnalysis.core.universe.Universe
New universe containing the original target atoms followed by the
inserted projectile copies.
See Also
--------
InsertPlanar : Slower per-molecule variant with full placement control.
SolvateCylinder : Equivalent for cylindrical regions.
"""
# Use no fewer than 20 atoms for solvation
SOLVATION_THRESHOLD = 20
if xmax is None:
xmax = TargetUniverse.dimensions[0]
if ymax is None:
ymax = TargetUniverse.dimensions[1]
if zmax is None:
zmax = TargetUniverse.dimensions[2]
if xmin is None:
xmin = 0
if ymin is None:
ymin = 0
if zmin is None:
zmin = 0
InsertionDomain = np.array([xmin, ymin, zmin, xmax, ymax, zmax])
for i in np.arange(3):
if InsertionDomain[i + 3] is None:
InsertionDomain[i + 3] = TargetUniverse.dimensions[i]
InsertionDomainSize = InsertionDomain[3:6] - InsertionDomain[0:3]
dimensionsTarget = TargetUniverse.dimensions.copy()
if density is not None:
n = np.floor(
density
* InsertionDomainSize[0]
* InsertionDomainSize[1]
* InsertionDomainSize[2]
)
nAtomsTarget = TargetUniverse.atoms.n_atoms
nAtomsProjectile = ProjectileUniverse.atoms.n_atoms
logger.info(f"Should solvate {n} Projectiles")
x = np.ceil((n / (solvate_factor * fudge_factor)) ** (1 / 3)).astype(int)
if x <= 1:
x = 1
logger.info(f"Solvation factor: {solvate_factor}")
logger.info(f"Best tiling is {x}x{x}x{x}.")
return _renumber_projectile_resids(
InsertPlanar(
TargetUniverse,
ProjectileUniverse,
n,
xmin,
ymin,
zmin,
xmax,
ymax,
zmax,
distance,
tries,
),
nAtomsTarget,
)
if n / (x**3) < SOLVATION_THRESHOLD and x > 2:
x -= 1
real_solvate_factor = n / (x**3)
logger.info(f"Solvation factor: {solvate_factor}")
logger.info(f"Best tiling is {x}x{x}x{x}.")
real_solvate_factor = np.ceil(real_solvate_factor * fudge_factor).astype(int)
logger.info("Real solvation factor is", real_solvate_factor)
logger.info(
"This results in a total of",
x**3 * (real_solvate_factor),
"projectiles in the solvate box",
)
solvate_box_dimensions = np.concatenate(
[InsertionDomainSize / x, dimensionsTarget[3:6]]
)
solvate_box = InsertPlanar(
empty(solvate_box_dimensions),
ProjectileUniverse,
real_solvate_factor,
distance=distance,
tries=tries * 1000,
)
# We tile the small box to make a big box that is big enough to contain
# the insertion domain
logger.info("Tiling solvate box...")
big_solvate_box = tile_universe(solvate_box, x, x, x)
# Shift the solvate box to the beginning of the insertion domain
big_solvate_box.atoms.translate(InsertionDomain[0:3])
logger.info("Inserting solvate box into target universe...")
nAtomsSolvate = big_solvate_box.atoms.n_atoms
logger.info("Target atoms:", nAtomsTarget)
logger.info("Projectile atoms:", nAtomsSolvate)
if nAtomsTarget == 0:
SolvatedUniverse = big_solvate_box
else:
SolvatedUniverse = mda.Merge(TargetUniverse.atoms, big_solvate_box.atoms)
SolvatedUniverse.dimensions = dimensionsTarget
target = SolvatedUniverse.atoms[0:nAtomsTarget]
projectile = SolvatedUniverse.atoms[-nAtomsSolvate:]
logger.info("Search for overlapping atoms...")
ns = mda.lib.NeighborSearch.AtomNeighborSearch(
projectile, SolvatedUniverse.dimensions
)
touching_atoms = ns.search(target, distance, level="R").atoms
if touching_atoms.n_atoms > 0:
# touching_atoms = touching_atoms.intersection(projectile).residues.atoms
# if touching_atoms.n_atoms / nAtomsProjectile:
logger.info(
"Removing touching projectiles:", touching_atoms.n_atoms / nAtomsProjectile
)
SolvatedUniverse = mda.Merge(SolvatedUniverse.atoms - touching_atoms)
SolvatedUniverse.dimensions = dimensionsTarget
logger.info("Resulting number of atoms:", SolvatedUniverse.atoms.n_atoms)
logger.info("Expected number of atoms:", n * nAtomsProjectile + nAtomsTarget)
missingProjectiles = int(
((n * nAtomsProjectile + nAtomsTarget) - SolvatedUniverse.atoms.n_atoms)
/ nAtomsProjectile
)
if density is not None:
logger.info(
f" {SolvatedUniverse.atoms.n_atoms - nAtomsTarget} projectiles inserted"
)
return _renumber_projectile_resids(SolvatedUniverse, nAtomsTarget)
if missingProjectiles > 0:
logger.info("Missing", missingProjectiles, "Projectiles.")
logger.info("Adjusting fudge factor and trying again.")
return SolvatePlanar(
TargetUniverse,
ProjectileUniverse,
n,
density,
xmin,
ymin,
zmin,
xmax,
ymax,
zmax,
distance,
solvate_factor,
fudge_factor + 10 * missingProjectiles / n,
tries,
)
if missingProjectiles < 0:
nonTargetAtoms = SolvatedUniverse.atoms[nAtomsTarget:]
logger.info("Too many projectiles inserted:", -missingProjectiles)
logger.info(nonTargetAtoms.n_atoms)
logger.info(nonTargetAtoms.residues.n_residues)
logger.info(np.unique(nonTargetAtoms.residues.resids).shape)
logger.info("Removing", -missingProjectiles, "randomly selected projectiles.")
ToBeRemoved = nonTargetAtoms.residues[
np.random.choice(
np.arange(len(nonTargetAtoms.residues)),
-missingProjectiles,
replace=False,
)
]
SolvatedUniverse = mda.Merge(SolvatedUniverse.atoms - ToBeRemoved.atoms)
SolvatedUniverse.dimensions = dimensionsTarget
logger.info("Final number of atoms:", SolvatedUniverse.atoms.n_atoms)
return _renumber_projectile_resids(SolvatedUniverse, nAtomsTarget)
logger.info("All projectiles inserted correctly")
return _renumber_projectile_resids(SolvatedUniverse, nAtomsTarget)
[docs]
def InsertPlanar(
TargetUniverse: mda.Universe,
ProjectileUniverse: mda.Universe,
n: int = 1,
xmin: int = 0,
ymin: int = 0,
zmin: int = 0,
xmax: Optional[float] = None,
ymax: Optional[float] = None,
zmax: Optional[float] = None,
distance: float = 1.25,
tries: int = 1000,
) -> mda.Universe:
"""Insert ``n`` copies of a projectile into a rectangular region.
Each projectile is placed at a random position and orientation inside the
axis-aligned box defined by ``(xmin, ymin, zmin)`` and
``(xmax, ymax, zmax)``. Up to ``tries`` placement attempts are made per
projectile; a :class:`RuntimeError` is raised if no overlap-free position
is found.
Parameters
----------
TargetUniverse : MDAnalysis.core.universe.Universe
Universe to insert into. May be empty; ``TargetUniverse.dimensions``
then defines the simulation cell of the returned universe.
ProjectileUniverse : MDAnalysis.core.universe.Universe
Molecule that is inserted repeatedly.
n : int, default 1
Number of projectile copies to insert.
xmin, ymin, zmin : float, default 0
Lower bounds of the insertion box, in Å.
xmax, ymax, zmax : float, optional
Upper bounds of the insertion box, in Å. Each defaults to the
corresponding component of ``TargetUniverse.dimensions``.
distance : float, default 1.25
Minimum allowed distance (Å) between the inserted projectile and any
existing atom in the target.
tries : int, default 1000
Maximum number of random placement attempts per projectile.
Returns
-------
MDAnalysis.core.universe.Universe
New universe containing the target atoms followed by the inserted
projectile copies.
Raises
------
RuntimeError
If no overlap-free position is found within ``tries`` attempts for a
given projectile.
See Also
--------
SolvatePlanar : Fast variant for many projectiles.
InsertCylinder, InsertSphere
"""
nAtomsTargetOriginal = TargetUniverse.atoms.n_atoms
InsertionDomain = [xmin, ymin, zmin, xmax, ymax, zmax]
for i in np.arange(3):
if InsertionDomain[i + 3] is None:
InsertionDomain[i + 3] = TargetUniverse.dimensions[i]
InsertionDomain = np.array(InsertionDomain)
nAtomsProjectile = ProjectileUniverse.atoms.n_atoms
dimensionsTarget = TargetUniverse.dimensions.copy()
ProjectileUniverse.atoms.translate(-ProjectileUniverse.atoms.center_of_geometry())
if TargetUniverse.atoms.n_atoms == 0:
TargetUniverse = ProjectileUniverse.copy()
TargetUniverse.dimensions = dimensionsTarget
TargetUniverse.atoms.translate(
pos_random(InsertionDomain) - ProjectileUniverse.atoms.center_of_geometry()
)
TargetUniverse.atoms.rotateby(*rot_random())
n -= 1
for _N in tqdm(np.arange(n)):
nAtomsTarget = TargetUniverse.atoms.n_atoms
TargetUniverse = mda.Merge(TargetUniverse.atoms, ProjectileUniverse.atoms)
TargetUniverse.dimensions = dimensionsTarget
target = TargetUniverse.atoms[0:nAtomsTarget]
projectile = TargetUniverse.atoms[-nAtomsProjectile:]
ns = mda.lib.NeighborSearch.AtomNeighborSearch(target, dimensionsTarget)
for _attempt in range(tries):
projectile.translate(
pos_random(InsertionDomain) - projectile.atoms.center_of_geometry()
)
projectile.rotateby(*rot_random())
if len(ns.search(projectile, distance)) == 0:
break
else:
raise RuntimeError(
"Error: No suitable position found,\
maybe you are trying to insert to many particles? Aborting."
)
return _renumber_projectile_resids(TargetUniverse, nAtomsTargetOriginal)
[docs]
def InsertCylinder(
TargetUniverse: mda.Universe,
ProjectileUniverse: mda.Universe,
n: int = 1,
pos: Optional[np.ndarray] = None,
radius: Optional[float] = None,
min: float = 0,
max: Optional[float] = None,
dim: int = 2,
distance: float = 1.25,
tries: int = 1000,
) -> mda.Universe:
"""Insert ``n`` copies of a projectile into a cylindrical region.
Each projectile is placed at a random position and orientation inside the
cylinder centred at ``pos`` with radius ``radius``, extending from ``min``
to ``max`` along axis ``dim``. Up to ``tries`` placement attempts are
made per projectile; a :class:`RuntimeError` is raised if no overlap-free
position is found.
Parameters
----------
TargetUniverse : MDAnalysis.core.universe.Universe
Universe to insert into. May be empty.
ProjectileUniverse : MDAnalysis.core.universe.Universe
Molecule that is inserted repeatedly.
n : int, default 1
Number of projectile copies to insert.
pos : array_like of shape (3,), optional
Centre of the cylinder, in Å. Defaults to the centre of geometry of
``TargetUniverse`` (or the centre of its box if the target is empty).
The component along ``dim`` is overridden by ``min``.
radius : float, optional
Cylinder radius, in Å. Defaults to half of the smallest box edge of
``TargetUniverse``.
min, max : float, optional
Lower and upper bound of the cylinder along axis ``dim``, in Å.
``max`` defaults to ``TargetUniverse.dimensions[dim]``.
dim : {0, 1, 2}, default 2
Index of the axis along which the cylinder extends (0 = x, 1 = y,
2 = z).
distance : float, default 1.25
Minimum allowed distance (Å) between the inserted projectile and any
existing atom in the target.
tries : int, default 1000
Maximum number of random placement attempts per projectile.
Returns
-------
MDAnalysis.core.universe.Universe
New universe containing the target atoms followed by the inserted
projectile copies.
Raises
------
RuntimeError
If no overlap-free position is found within ``tries`` attempts for a
given projectile.
See Also
--------
SolvateCylinder : Fast variant for many projectiles.
InsertPlanar, InsertSphere
"""
nAtomsTargetOriginal = TargetUniverse.atoms.n_atoms
if max is None:
max = TargetUniverse.dimensions[dim]
nAtomsProjectile = ProjectileUniverse.atoms.n_atoms
dimensionsTarget = TargetUniverse.dimensions.copy()
if pos is None:
if TargetUniverse.atoms.n_atoms == 0:
pos = dimensionsTarget / 2
else:
pos = TargetUniverse.atoms.center_of_geometry()
pos[dim] = min
if radius is None:
radius = np.min(dimensionsTarget) / 2
ProjectileUniverse.atoms.translate(-ProjectileUniverse.atoms.center_of_geometry())
for _N in tqdm(np.arange(n)):
nAtomsTarget = TargetUniverse.atoms.n_atoms
TargetUniverse = mda.Merge(TargetUniverse.atoms, ProjectileUniverse.atoms)
TargetUniverse.dimensions = dimensionsTarget.copy()
target = TargetUniverse.atoms[0:nAtomsTarget]
projectile = TargetUniverse.atoms[-nAtomsProjectile:]
ns = mda.lib.NeighborSearch.AtomNeighborSearch(target)
# Generate coordinates and check for overlap
for _attempt in range(tries):
projectile.rotateby(*rot_random())
r = radius * np.sqrt(np.random.rand())
phi, z = np.random.rand(2) * [2 * np.pi, (max - min)]
newcoord = np.roll([r * np.cos(phi), r * np.sin(phi), z], dim - 2) + pos
projectile.translate(newcoord - projectile.atoms.center_of_geometry())
if len(ns.search(projectile, distance)) == 0:
break
else:
raise RuntimeError(
"Error: No suitable position found,\
maybe you are trying to insert too many particles? Aborting."
)
return _renumber_projectile_resids(TargetUniverse, nAtomsTargetOriginal)
[docs]
def InsertSphere(
TargetUniverse: mda.Universe,
ProjectileUniverse: mda.Universe,
n: int = 1,
pos: Optional[np.ndarray] = None,
radius: Optional[float] = None,
distance: float = 1.25,
tries: int = 1000,
) -> mda.Universe:
"""Insert ``n`` copies of a projectile into a spherical region.
Each projectile is placed at a uniformly random position inside the
sphere centred at ``pos`` with radius ``radius`` and a random
orientation. Up to ``tries`` placement attempts are made per projectile;
a :class:`RuntimeError` is raised if no overlap-free position is found.
Parameters
----------
TargetUniverse : MDAnalysis.core.universe.Universe
Universe to insert into. May be empty.
ProjectileUniverse : MDAnalysis.core.universe.Universe
Molecule that is inserted repeatedly.
n : int, default 1
Number of projectile copies to insert.
pos : array_like of shape (3,), optional
Centre of the sphere, in Å. Defaults to the centre of geometry of
``TargetUniverse`` (or the centre of its box if the target is empty).
radius : float, optional
Sphere radius, in Å. Defaults to half of the smallest box edge of
``TargetUniverse``.
distance : float, default 1.25
Minimum allowed distance (Å) between the inserted projectile and any
existing atom in the target.
tries : int, default 1000
Maximum number of random placement attempts per projectile.
Returns
-------
MDAnalysis.core.universe.Universe
New universe containing the target atoms followed by the inserted
projectile copies.
Raises
------
RuntimeError
If no overlap-free position is found within ``tries`` attempts for a
given projectile.
See Also
--------
InsertPlanar, InsertCylinder
"""
def rand_spherical(radius: float = 1.0) -> np.ndarray:
u = np.random.rand()
v = np.random.rand()
theta = u * 2.0 * np.pi
phi = np.arccos(2.0 * v - 1.0)
r = radius * np.power(np.random.rand(), 1 / 3)
sinTheta = np.sin(theta)
cosTheta = np.cos(theta)
sinPhi = np.sin(phi)
cosPhi = np.cos(phi)
x = r * sinPhi * cosTheta
y = r * sinPhi * sinTheta
z = r * cosPhi
return np.array([x, y, z])
nAtomsTargetOriginal = TargetUniverse.atoms.n_atoms
nAtomsProjectile = ProjectileUniverse.atoms.n_atoms
dimensionsTarget = TargetUniverse.dimensions.copy()
if pos is None:
if TargetUniverse.atoms.n_atoms == 0:
pos = dimensionsTarget[:3] / 2
else:
pos = TargetUniverse.atoms.center_of_geometry()
if radius is None:
radius = np.min(dimensionsTarget) / 2
ProjectileUniverse.atoms.translate(-ProjectileUniverse.atoms.center_of_geometry())
if TargetUniverse.atoms.n_atoms == 0:
TargetUniverse = ProjectileUniverse.copy()
TargetUniverse.dimensions = dimensionsTarget
TargetUniverse.atoms.translate(
pos + rand_spherical(radius) - TargetUniverse.atoms.center_of_geometry()
)
TargetUniverse.atoms.rotateby(*rot_random())
n -= 1
for _N in tqdm(np.arange(n)):
nAtomsTarget = TargetUniverse.atoms.n_atoms
TargetUniverse = mda.Merge(TargetUniverse.atoms, ProjectileUniverse.atoms)
TargetUniverse.dimensions = dimensionsTarget.copy()
target = TargetUniverse.atoms[0:nAtomsTarget]
projectile = TargetUniverse.atoms[-nAtomsProjectile:]
ns = mda.lib.NeighborSearch.AtomNeighborSearch(target)
# Generate coordinates and check for overlap
for _attempt in range(tries):
projectile.rotateby(*rot_random())
newcoord = rand_spherical(radius) + pos
projectile.translate(newcoord - projectile.atoms.center_of_geometry())
if len(ns.search(projectile, distance)) == 0:
break
else:
raise RuntimeError(
"Error: No suitable position found, \
maybe you are trying to insert to many particles? Aborting."
)
return _renumber_projectile_resids(TargetUniverse, nAtomsTargetOriginal)