Skip to content

Commit

Permalink
Merge branch 'python' into fix-3119
Browse files Browse the repository at this point in the history
  • Loading branch information
@kodiakhq
kodiakhq[bot] authored Feb 4, 2020
2 parents 7f7dda5 + 780c7fd commit 4047fdb
Show file tree
Hide file tree
Showing 11 changed files with 624 additions and 116 deletions.
1 change: 1 addition & 0 deletions requirements.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,3 +27,4 @@ pycodestyle==2.3.1
pylint>=2.2.2
astroid>=2.1.0
isort>=4.3.4
setuptools>=20.7.0
110 changes: 58 additions & 52 deletions samples/rigid_body.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,42 +19,35 @@
``VIRTUAL_SITES_RELATIVE`` feature.
"""

import espressomd
required_features = ["VIRTUAL_SITES_RELATIVE", "MASS", "ROTATIONAL_INERTIA"]
espressomd.assert_features(required_features)
from espressomd.virtual_sites import VirtualSitesRelative
import enum
import math

import numpy as np

import espressomd
required_features = ["VIRTUAL_SITES_RELATIVE", "MASS", "ROTATIONAL_INERTIA"]
espressomd.assert_features(required_features)
import espressomd.virtual_sites
import espressomd.rotation

box_l = 100
system = espressomd.System(box_l=[box_l, box_l, box_l])
system.set_random_state_PRNG()
system.seed = system.cell_system.get_state()['n_nodes'] * [1234]

system = espressomd.System(box_l=[10.0] * 3)
system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative(
have_velocity=True)
system.time_step = 0.01
skin = 10.0
system.cell_system.skin = skin
system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=42)

# Particle types
type_centre = 0
type_A = 1

system.thermostat.set_langevin(kT=1.0, gamma=20.0, seed=42)

#############################################################
print("** Placing particles")
#############################################################

class ParticleTypes(enum.IntEnum):
CENTER = enum.auto()
BRANCH = enum.auto()

# start at p_id=1, and reserve p_id=0 for the centre bead.
p_id = 1

branch_len = 5
x = y = z = box_l / 2
center = 0.5 * system.box_l

# place six branches, pointing +/-x +/-y and +/-z
# note that we do not make the particles virtual at this point.
# Place six branches, pointing +/-x +/-y and +/-z.
# Note that we do not make the particles virtual at this point.
# The script uses center of mass an moment of inertia analysis routines
# to obtain the position and inertia moments of the central particle.
# Once a particle is made virtual, it will no longer contribute to
Expand All @@ -63,47 +56,60 @@

for direction in np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]):
for n in range(branch_len):
system.part.add(pos=[x, y, z] + (n + 1) * direction,
type=type_A)
system.part.add(pos=[x, y, z] - (n + 1) * direction,
type=type_A, virtual=True)
system.part.add(pos=center + (n + 1) * direction,
type=ParticleTypes.BRANCH.value)
system.part.add(pos=center - (n + 1) * direction,
type=ParticleTypes.BRANCH.value)

system.virtual_sites = VirtualSitesRelative(have_velocity=True)

# here we calculate the center of mass position (com) and the moment of
# inertia (momI) of the object
com = system.analysis.center_of_mass(p_type=type_A)
print("center of mass is:", com)
center_of_mass = system.analysis.center_of_mass(
p_type=ParticleTypes.BRANCH.value)
print("Center of mass:", center_of_mass)

# if using multiple nodes, we need to change min_global_cut to the largest
# separation
if system.cell_system.get_state()['n_nodes'] > 1:
max_dist = 0
for p in system.part:
dist = np.sum((p.pos - com)**2)
if dist > max_dist:
max_dist = dist
max_dist = np.sqrt(max_dist)
print("max dist is {}".format(max_dist))
system.min_global_cut = max_dist

mat_I = system.analysis.moment_of_inertia_matrix(p_type=type_A)
max_inter = np.max(np.linalg.norm(system.part[:].pos - center_of_mass, axis=1))
system.min_global_cut = max_inter


principal_moments, principal_axes = espressomd.rotation.diagonalized_inertia_tensor(
system.part[:].pos, system.part[:].mass)
# in this simple case, the cluster has principal axes aligned with the box
momI = [mat_I[0, 0], mat_I[1, 1], mat_I[2, 2]]
print("moment of inertia is", momI)
print("Principal moments: {}, principal axes tensor: {}".format(
principal_moments, principal_axes))

# if we rotate the arms, we have to make sure that we set the quaternion of the
# center particle accordingly while setting the principal moments of inertia
AXIS = np.array([1., 0., 0.])
ANGLE = np.pi / 4.0


def rotate_vector(vector, axis, angle):
return axis * np.dot(axis, vector) + math.cos(angle) * np.cross(
np.cross(axis, vector), axis) + math.sin(angle) * np.cross(axis, vector)


for p in system.part:
p.pos = rotate_vector(p.pos - center_of_mass, AXIS, ANGLE) + center_of_mass


principal_moments, principal_axes = espressomd.rotation.diagonalized_inertia_tensor(
system.part[:].pos, system.part[:].mass)
# after rotating the whole object the principal axes changed
print("After rotating: principal moments: {}, principal axes tensor: {}".format(
principal_moments, principal_axes))

# place center bead
p_center = system.part.add(
pos=com, mass=branch_len * 6 + 1, rinertia=momI,
rotation=[1, 1, 1], type=type_centre)
pos=center_of_mass, mass=branch_len * 6 + 1, rinertia=principal_moments,
rotation=[1, 1, 1], type=ParticleTypes.CENTER.value, quat=espressomd.rotation.matrix_to_quat(principal_axes))

# Relate the particles that make up the rigid body to the central particle.
# This will also mark them as `virtual = True`
for p in system.part:
if p != p_center:
p.vs_auto_relate_to(p_center.id)
for p in system.part.select(type=ParticleTypes.BRANCH.value):
p.vs_auto_relate_to(p_center.id)

for frame in range(200):
system.integrator.run(100)

print("**Simulation finished")
print("Simulation finished")
116 changes: 72 additions & 44 deletions src/core/bonded_interactions/bonded_interaction_data.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,59 +19,86 @@
#include "bonded_interaction_data.hpp"
#include "communication.hpp"

#include <boost/range/numeric.hpp>
#include <utils/constants.hpp>

std::vector<Bonded_ia_parameters> bonded_ia_params;

auto cutoff(int type, Bond_parameters const &bp) {
switch (type) {
case BONDED_IA_NONE:
return -1.;
case BONDED_IA_FENE:
return bp.fene.cutoff();
case BONDED_IA_HARMONIC:
return bp.harmonic.cutoff();
case BONDED_IA_HARMONIC_DUMBBELL:
return bp.harmonic_dumbbell.cutoff();
case BONDED_IA_QUARTIC:
return bp.quartic.cutoff();
case BONDED_IA_BONDED_COULOMB:
return bp.bonded_coulomb.cutoff();
case BONDED_IA_BONDED_COULOMB_SR:
return bp.bonded_coulomb_sr.cutoff();
case BONDED_IA_DIHEDRAL:
return bp.dihedral.cutoff();
case BONDED_IA_TABULATED_DISTANCE:
case BONDED_IA_TABULATED_ANGLE:
case BONDED_IA_TABULATED_DIHEDRAL:
return bp.tab.cutoff();
case BONDED_IA_SUBT_LJ:
return bp.subt_lj.cutoff();
case BONDED_IA_RIGID_BOND:
return bp.rigid_bond.cutoff();
case BONDED_IA_VIRTUAL_BOND:
return bp.virt.cutoff();
case BONDED_IA_ANGLE_HARMONIC:
return bp.angle_harmonic.cutoff();
case BONDED_IA_ANGLE_COSINE:
return bp.angle_cosine.cutoff();
case BONDED_IA_ANGLE_COSSQUARE:
return bp.angle_cossquare.cutoff();
case BONDED_IA_OIF_LOCAL_FORCES:
return bp.oif_local_forces.cutoff();
case BONDED_IA_OIF_GLOBAL_FORCES:
return bp.oif_global_forces.cutoff();
case BONDED_IA_IBM_TRIEL:
return bp.ibm_triel.cutoff();
case BONDED_IA_IBM_VOLUME_CONSERVATION:
return bp.ibmVolConsParameters.cutoff();
case BONDED_IA_IBM_TRIBEND:
return bp.ibm_tribend.cutoff();
case BONDED_IA_UMBRELLA:
return bp.umbrella.cutoff();
case BONDED_IA_THERMALIZED_DIST:
return bp.thermalized_bond.cutoff();
default:
throw std::runtime_error("Unknown bond type.");
}
}

double recalc_maximal_cutoff_bonded() {
auto max_cut_bonded = -1.;
auto const max_cut_bonded =
boost::accumulate(bonded_ia_params, -1.,
[](auto max_cut, Bonded_ia_parameters const &bond) {
return std::max(max_cut, cutoff(bond.type, bond.p));
});

for (auto const &bonded_ia_param : bonded_ia_params) {
switch (bonded_ia_param.type) {
case BONDED_IA_FENE:
max_cut_bonded =
std::max(max_cut_bonded,
bonded_ia_param.p.fene.r0 + bonded_ia_param.p.fene.drmax);
break;
case BONDED_IA_HARMONIC:
max_cut_bonded =
std::max(max_cut_bonded, bonded_ia_param.p.harmonic.r_cut);
break;
case BONDED_IA_THERMALIZED_DIST:
max_cut_bonded =
std::max(max_cut_bonded, bonded_ia_param.p.thermalized_bond.r_cut);
break;
case BONDED_IA_RIGID_BOND:
max_cut_bonded =
std::max(max_cut_bonded, std::sqrt(bonded_ia_param.p.rigid_bond.d2));
break;
case BONDED_IA_TABULATED_DISTANCE:
max_cut_bonded =
std::max(max_cut_bonded, bonded_ia_param.p.tab.pot->cutoff());
break;
case BONDED_IA_IBM_TRIEL:
max_cut_bonded =
std::max(max_cut_bonded, bonded_ia_param.p.ibm_triel.maxDist);
break;
default:
break;
}
}
/* Check if there are dihedrals */
auto const any_dihedrals = std::any_of(
bonded_ia_params.begin(), bonded_ia_params.end(), [](auto const &bond) {
switch (bond.type) {
case BONDED_IA_DIHEDRAL:
case BONDED_IA_TABULATED_DIHEDRAL:
return true;
default:
return false;
}
});

/* dihedrals: the central particle is indirectly connected to the fourth
* particle via the third particle, so we have to double the cutoff */
for (auto const &bonded_ia_param : bonded_ia_params) {
switch (bonded_ia_param.type) {
case BONDED_IA_DIHEDRAL:
case BONDED_IA_TABULATED_DIHEDRAL:
max_cut_bonded *= 2;
break;
default:
break;
}
}

return max_cut_bonded;
return (any_dihedrals) ? 2 * max_cut_bonded : max_cut_bonded;
}

void make_bond_type_exist(int type) {
Expand All @@ -95,6 +122,7 @@ int virtual_set_params(int bond_type) {

bonded_ia_params[bond_type].type = BONDED_IA_VIRTUAL_BOND;
bonded_ia_params[bond_type].num = 1;
bonded_ia_params[bond_type].p.virt = VirtualBond_Parameters{};

/* broadcast interaction parameters */
mpi_bcast_ia_params(bond_type, -1);
Expand Down
Loading

0 comments on commit 4047fdb

Please sign in to comment.