# Copyright 2021-2023 @ Shenzhen Bay Laboratory &
# Peking University &
# Huawei Technologies Co., Ltd
#
# This code is a part of MindSPONGE:
# MindSpore Simulation Package tOwards Next Generation molecular modelling.
#
# MindSPONGE is open-source software based on the AI-framework:
# MindSpore (https://www.mindspore.cn/)
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ============================================================================
"""
Common operations
"""
from inspect import signature
import numpy as np
import mindspore as ms
from mindspore import numpy as msnp
from mindspore.ops import functional as F
from mindspore import ops, nn
from mindspore import Tensor
from mindspore.nn import Cell
from . import functions as func
from .functions import get_integer
from .units import Units, GLOBAL_UNITS
__all__ = [
'GetVector',
'GetDistance',
'VelocityGenerator',
'GetDistanceShift',
'GetShiftGrad',
]
[docs]class GetVector(Cell):
r"""The class to get vector with or without PBC box
Args:
use_pbc (bool): Whether to calculate vector under periodic boundary condition.
If ``None`` is given, it will determine whether to use periodic boundary
conditions based on whether the ``pbc_box`` is provided.
Default: ``None`` .
Examples:
>>> import mindspore as ms
>>> import numpy as np
>>> from sponge.function import GetVector
>>> from mindspore import Tensor
>>> crd = Tensor(np.random.random((4, 3)), ms.float32)
>>> pbc_box = Tensor([[0.5, 0.5, 0.5]], ms.float32)
>>> gd = GetVector(use_pbc=True)
>>> gd(crd[0], crd[1], pbc_box)
Tensor(shape=[1, 3], dtype=Float32, value=
[[ 5.14909625e-02, -4.62748706e-02, -1.20763242e-01]])
>>> gd = GetVector(use_pbc=False)
>>> gd(crd[0], crd[1])
Tensor(shape=[3], dtype=Float32, value= [ 5.14909625e-02, 4.53725129e-01, -1.20763242e-01])
"""
def __init__(self, use_pbc: bool = None):
super().__init__()
self._use_pbc = use_pbc
if use_pbc is None:
self.calc_vector = self.calc_vector_default
else:
if use_pbc:
self.calc_vector = self.calc_vector_pbc
else:
self.calc_vector = self.calc_vector_nopbc
@property
def use_pbc(self) -> bool:
"""whether to use periodic boundary condition
Returns:
bool, whether to use periodic boundary condition
"""
return self._use_pbc
@use_pbc.setter
def use_pbc(self, use_pbc_: bool):
"""set whether to use periodic boundary condition"""
self.set_pbc(use_pbc_)
[docs] def set_pbc(self, use_pbc: bool):
"""
set whether to use periodic boundary condition.
Args:
use_pbc (bool): Whether to calculate vector under periodic boundary condition.
Default: ``None``.
"""
self._use_pbc = use_pbc
if use_pbc is None:
self.calc_vector = self.calc_vector_default
else:
if not isinstance(use_pbc, bool):
raise TypeError(f'The type of use_pbc must be bool or None but got: {type(use_pbc)}')
if use_pbc:
self.calc_vector = self.calc_vector_pbc
else:
self.calc_vector = self.calc_vector_nopbc
return self
[docs] def calc_vector_default(self, initial: Tensor, terminal: Tensor, pbc_box: Tensor = None) -> Tensor:
"""
get vector.
Args:
initial (Tensor): Tensor of shape :math:`(B, ..., D)` .
B means batchsize, i.e. number of walkers in simulation.
D means spatial dimension of the simulation system. Usually is 3.
Data type is float.
Coordinate of initial point
terminal (Tensor): Tensor of shape :math:`(B, ..., D)` . Data type is float.
Coordinate of terminal point
pbc_box (Tensor): Tensor of shape :math:`(B, D)` . Data type is float.
Default: ``None``.
"""
return func.calc_vector(initial, terminal, pbc_box)
[docs] def calc_vector_pbc(self, initial: Tensor, terminal: Tensor, pbc_box: Tensor = None) -> Tensor:
"""
get vector with perodic bundary condition.
Args:
initial (Tensor): Tensor of shape :math:`(B, ..., D)` .
B means batchsize, i.e. number of walkers in simulation.
D means spatial dimension of the simulation system. Usually is 3.
Data type is float.
Coordinate of initial point
terminal (Tensor): Tensor of shape :math:`(B, ..., D)` . Data type is float.
Coordinate of terminal point
pbc_box (Tensor): Tensor of shape :math:`(B, D)` . Data type is float.
Default: ``None``.
"""
return func.calc_vector_pbc(initial, terminal, pbc_box)
[docs] def calc_vector_nopbc(self, initial: Tensor, terminal: Tensor, pbc_box: Tensor = None) -> Tensor:
"""
get vector without perodic bundary condition.
Args:
initial (Tensor): Tensor of shape :math:`(B, ..., D)` .
B means batchsize, i.e. number of walkers in simulation.
D means spatial dimension of the simulation system. Usually is 3.
Data type is float.
Coordinate of initial point
terminal (Tensor): Tensor of shape :math:`(B, ..., D)` . Data type is float.
Coordinate of terminal point
pbc_box (Tensor): Tensor of shape :math:`(B, D)` . Data type is float.
Default: ``None``.
"""
#pylint: disable=unused-argument
return terminal - initial
def construct(self, initial: Tensor, terminal: Tensor, pbc_box: Tensor = None):
r"""Compute vector from initial point to terminal point.
Args:
initial (Tensor): Tensor of shape :math:`(B, ..., D)` .
B means batchsize, i.e. number of walkers in simulation.
D means spatial dimension of the simulation system. Usually is 3.
Data type is float.
Coordinate of initial point
terminal (Tensor): Tensor of shape :math:`(B, ..., D)` . Data type is float.
Coordinate of terminal point
pbc_box (Tensor): Tensor of shape :math:`(B, D)` . Data type is float.
Default: ``None``.
Returns:
vector (Tensor): Tensor of shape :math:`(B, ..., D)` . Data type is float.
"""
return self.calc_vector(initial, terminal, pbc_box)
[docs]class GetDistance(GetVector):
r"""The class to calculate distance with or without PBC box
Args:
use_pbc (bool): Whether to calculate distance under periodic boundary condition.
If this is "None", it will determine whether to calculate the distance under
periodic boundary condition based on whether the pbc_box is given.
Default: ``None``.
keepdims (bool): Whether to keep the last dimension of the output Tensor of distance after norm.
If this is "True", the last dimension of the output Tensor will be 1.
Default: ``False``.
axis (int): The axis of the space dimension of the coordinate. Default: ``-1`` .
Examples:
>>> import mindspore as ms
>>> import numpy as np
>>> from sponge.function import GetDistance
>>> from mindspore import Tensor
>>> crd = Tensor(np.random.random((4, 3)), ms.float32)
>>> pbc_box = Tensor([[0.5, 0.5, 0.5]], ms.float32)
>>> gd = GetDistance(use_pbc=True, keepdims=False)
>>> gd(crd[0], crd[1], pbc_box)
Tensor(shape=[1], dtype=Float32, value= [ 1.39199302e-01])
>>> gd = GetDistance(use_pbc=False, keepdims=False)
>>> gd(crd[0], crd[1])
Tensor(shape=[], dtype=Float32, value= 0.472336)
"""
def __init__(self,
use_pbc: bool = None,
keepdims: bool = False,
axis: int = -1,
):
super().__init__(use_pbc=use_pbc)
self.axis = get_integer(axis)
self.keepdims = keepdims
self.norm = None
# MindSpore < 2.0.0-rc1
if 'ord' not in signature(ops.norm).parameters.keys():
self.norm = nn.Norm(self.axis, self.keepdims)
def construct(self, initial: Tensor, terminal: Tensor, pbc_box: Tensor = None):
r"""Compute the distance from initial point to terminal point.
Args:
initial (Tensor): Tensor of shape (B, ..., D). Data type is float.
Coordinate of initial point
B means batchsize, i.e. number of walkers in simulation.
D means spatial dimension of the simulation system. Usually is 3.
terminal (Tensor): Tensor of shape (B, ..., D). Data type is float.
Coordinate of terminal point
pbc_box (Tensor): Tensor of shape (B, D). Data type is float.
Default: ``None``.
Returns:
distance (Tensor): Tensor of shape (B, ...). Data type is float.
"""
if self._use_pbc:
vector = self.calc_vector(initial, terminal, pbc_box)
else:
vector = self.calc_vector(initial, terminal)
if self.norm is None:
return ops.norm(vector, None, self.axis, self.keepdims)
return self.norm(vector)
[docs]class VelocityGenerator(Cell):
r"""A class to generate velocities for atoms in system according to temperature
Args:
temperature (float): Temperature
remove_translation (bool): Whether to calculate distance under periodic boundary condition.
Default: ``True``
seed (int): Random seed for standard normal. Default: ``0``
seed2 (int): Random seed2 for standard normal. Default: ``0``
length_unit (str): Length unit. Default: ``None``
energy_unit (str): energy unit. Default: ``None``
Examples:
>>> from sponge import UpdaterMD
>>> from sponge.function import VelocityGenerator
>>> vgen = VelocityGenerator(300)
>>> velocity = vgen(system.shape, system.atom_mass)
>>> opt = UpdaterMD(system=system,
... time_step=1e-3,
... velocity=velocity,
... integrator='velocity_verlet',
... temperature=300,
... thermostat='langevin')
"""
#pylint: disable=invalid-name
def __init__(self,
temperature: float = 300,
remove_translation: bool = True,
seed: int = 0,
seed2: int = 0,
length_unit: str = None,
energy_unit: str = None,
):
super().__init__()
if length_unit is None and energy_unit is None:
self.units = GLOBAL_UNITS
else:
self.units = Units(length_unit, energy_unit)
self.temperature = Tensor(temperature, ms.float32).reshape(-1, 1, 1)
self.standard_normal = ops.StandardNormal(seed, seed2)
self.kb = Tensor(self.units.boltzmann, ms.float32)
self.kbT = self.kb * self.temperature
self.sigma = F.sqrt(self.kbT)
self.kinectics_unit_scale = self.units.kinetic_ref
self.remove_translation = remove_translation
self.identity = ops.Identity()
self.multi_temp = False
[docs] def set_temperature(self, temperature: float):
"""
set temperature.
Args:
temperature (float): temperature
"""
self.temperature = Tensor(temperature, ms.float32).reshape(-1, 1, 1)
self.multi_temp = False
if self.temperature is not None and self.temperature.size > 1:
self.multi_temp = False
return self
def construct(self, shape: tuple, atom_mass: Tensor, mask: Tensor = None):
r"""Randomly generate velocities for atoms in system.
Args:
shape (tuple): Shape of velocity
atom_mass (Tensor): Tensor of shape (B, A). Data type is float.
Atom mass in system.
B means batchsize, i.e. number of walkers in simulation.
A means number of atoms
mask (Tensor): Tensor of shape (B, A). Data type is bool.
Mask for atoms. Default: ``None``.
Returns:
velocity (Tensor): Tensor of shape (B, A, D). Data type is float.
D means spatial dimension of the simulation system. Usually is 3.
"""
# (B,A,1)
atom_mass = F.expand_dims(self.identity(atom_mass), -1)
inv_mass = msnp.reciprocal(atom_mass)
velocity_scale = self.sigma * \
msnp.sqrt(inv_mass / self.kinectics_unit_scale)
if mask is not None:
velocity_scale = msnp.where(
F.expand_dims(mask, -1), velocity_scale, 0)
velocity = self.standard_normal(shape) * velocity_scale
if self.remove_translation:
# (B,A,D) * (1,A,1)
momentum = atom_mass * velocity
# (1,1,1) or (B,1,1) <- (1,A,1) or (B,A,1)
dp = func.keepdims_mean(momentum, -2)
if mask is not None:
sp = func.keepdims_sum(momentum, -2)
n = func.keepdims_sum(F.cast(mask, ms.int32), -2)
dp = sp / n
# (B,A,D) - (B,1,D) = (B,A,D)
momentum -= dp
velocity = momentum * inv_mass
return velocity
[docs]class GetDistanceShift(Cell):
r"""Module for calculating B matrix whose dimensions are: C.
Args:
bonds (Tensor): Tensor of shape (C, 2). Data type is int.
Bonds need to be constraint.
num_atoms (int): Number of atoms in system.
num_walkers (int): Number of multiple walkers.
use_pbc (bool): Whether to use periodic boundary condition.
"""
def __init__(self,
bonds: Tensor,
num_atoms: int,
num_walkers: int = 1,
use_pbc: bool = None
):
super().__init__(auto_prefix=False)
# (C,2)
self.bonds = bonds
# (B,C,A)
shape = (num_walkers, bonds.shape[-2], num_atoms)
# (1,C,1)
bond0 = self.bonds[..., 0].reshape(1, -1, 1).asnumpy()
# (B,C,A) <- (B,A,1)
mask0 = np.zeros(shape)
np.put_along_axis(mask0, bond0, 1, axis=-1)
# (B,C,A,1)
self.mask0 = F.expand_dims(Tensor(mask0, ms.int32), -1)
# (1,C,1)
bond1 = self.bonds[..., 1].reshape(1, -1, 1).asnumpy()
# (B,C,A) <- (B,A,1)
mask1 = np.zeros(shape)
np.put_along_axis(mask1, bond1, 1, axis=-1)
# (B,C,A,1)
self.mask1 = F.expand_dims(Tensor(mask1, ms.int32), -1)
self.get_distance = GetDistance(use_pbc)
def construct(self, coordinate_new: Tensor, coordinate_old: Tensor, pbc_box: Tensor = None):
"""Module for calculating B matrix whose dimensions are: C.
Args:
coordinate_new (Tensor): Tensor of shape (B,A,D). Data type is float.
The new coordinates of the system.
coordinate_old (Tensor): Tensor of shape (B,A,D). Data type is float.
The old coordinates of the system.
pbc_box (Tensor): Tensor of shape (B,D). Data type is float.
Tensor of PBC box
Returns:
shift (Tensor): Tensor of shape (B,A,D). Data type is float.
"""
# (B,C,A,D) = (B,C,A,1) * (B,1,A,D)
pos_new_0 = F.reduce_sum(self.mask0 * coordinate_new, -2)
pos_new_1 = F.reduce_sum(self.mask1 * coordinate_new, -2)
# (B,C,A)
dis_new = self.get_distance(pos_new_0, pos_new_1, pbc_box)
# (B,C,A,D) = (B,C,A,1) * (B,1,A,D)
pos_old_0 = F.reduce_sum(self.mask0 * coordinate_old, -2)
pos_old_1 = F.reduce_sum(self.mask1 * coordinate_old, -2)
dis_old = self.get_distance(pos_old_0, pos_old_1, pbc_box)
# (B,C,A)
return dis_new - dis_old
[docs]class GetShiftGrad(Cell):
"""Module for calculating the differentiation of B matrix whose dimensions are: K*N*D.
Args:
num_atoms (int): Number of atoms in system.
bonds (Tensor): Tensor of shape :math:`(C, 2)` . Data type is int.
Bonds need to be constraint.
num_walkers (int): Number of multiple walkers.
dimension (int): Dimension.
use_pbc (bool): Whether to use periodic boundary condition.
"""
def __init__(self,
num_atoms: int,
bonds: Tensor,
num_walkers: int = 1,
dimension: int = 3,
use_pbc: bool = None
):
super().__init__(auto_prefix=False)
# (B,C,A,D)
shape = (num_walkers, bonds.shape[-2], num_atoms, dimension)
self.broadcast = ops.BroadcastTo(shape)
self.net = GetDistanceShift(
bonds=bonds,
num_atoms=num_atoms,
num_walkers=num_walkers,
use_pbc=use_pbc
)
self.grad = ops.GradOperation()
self.zero_shift = ops.Zeros()((num_walkers, num_atoms - 1, num_atoms, dimension), ms.float32)
def construct(self, coordinate_new: Tensor, coordinate_old: Tensor, pbc_box: Tensor = None):
"""Module for calculating the differentiation of B matrix whose dimensions are: K*N*D.
Args:
coordinate_new (Tensor): Tensor of shape :math:`(B,A,D)` . Data type is float.
The new coordinates of the system.
coordinate_old (Tensor): Tensor of shape :math:`(B,A,D)` . Data type is float.
The old coordinates of the system.
pbc_box (Tensor): Tensor of shape :math:`(B,D)` . Data type is float.
Tensor of PBC box
Returns:
shift (Tensor): Tensor of shape :math:`(B,A,D)` . Data type is float.
"""
# (B,C,A,D)
coordinate_new = self.broadcast(coordinate_new[:, None, :, :])
coordinate_old = self.broadcast(coordinate_old[:, None, :, :])
shift_grad = self.grad(self.net)(coordinate_new, coordinate_old, pbc_box)
if msnp.isnan(shift_grad.sum()):
shift_grad = self.zero_shift
return shift_grad