# Copyright (C) 2010-2018 The ESPResSo project
#
# This file is part of ESPResSo.
#
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from __future__ import print_function
import espressomd.interactions
from espressomd import has_features
# Dict with drude type infos
drude_dict = {}
# Lists with unique drude and core types
core_type_list = []
drude_type_list = []
# Get core id from drude id
core_id_from_drude_id = {}
# Drude IDs
drude_id_list = []
[docs]def add_drude_particle_to_core(system, harmonic_bond, thermalized_bond,
p_core, id_drude, type_drude, alpha,
mass_drude, coulomb_prefactor,
thole_damping=2.6, verbose=False):
"""
Adds a drude particle with specified id, type, and mass to the system.
Checks if different Drude particles have different types.
Collects types/charges/polarizations/Thole factors for intramol. core-Drude short-range exclusion and Thole interaction.
Attributes
----------
system : Instance of :attr:`espressomd.System`
harmonic_bond: This method adds this harmonic bond to between drude particle and core
thermalized_bond: This method adds this thermalizerd_bond to between drude particle and core
p_core: The existing core particle
id_drude: :obj:`int`
This method creates the drude particle and assigns this id.
type_drude: :obj:`int`
The type of the newly created drude particle
alpha : :obj:`float`
The polarizability in units of inverse volume. Related to the charge of the drude particle.
mass_drude : :obj:`float`
The mass of the newly created drude particle
coulomb_prefactor : :obj:`float`
Required to calculate the charge of the drude particle.
thole_damping : :obj:`float`
Thole damping factor of the drude pair. Comes to effect if add_all_thole() method is used.
verbose : :obj:`bool`
Turns on verbosity.
"""
k = harmonic_bond.params["k"]
q_drude = -1.0 * pow(k * alpha / coulomb_prefactor, 0.5)
if has_features("PARTICLE_ANISOTROPY"):
gamma_off = [0.0, 0.0, 0.0]
else:
gamma_off = 0.0
system.part.add(id=id_drude, pos=p_core.pos, type=type_drude,
q=q_drude, mass=mass_drude, temp=0, gamma=gamma_off)
if verbose:
print(
"Adding to core", p_core.id, "drude id", id_drude, " pol", alpha,
" core charge", p_core.q, "->", p_core.q - q_drude, " drude charge", q_drude)
p_core.q -= q_drude
p_core.mass -= mass_drude
p_core.add_bond((harmonic_bond, id_drude))
p_core.add_bond((thermalized_bond, id_drude))
p_core.temp = 0.0
p_core.gamma = gamma_off
if type_drude in drude_dict and not (drude_dict[type_drude]["q"] == q_drude and drude_dict[type_drude]["thole_damping"] == thole_damping):
raise Exception(
"Drude particles with different drude charges have to have different types for THOLE")
core_id_from_drude_id[id_drude] = p_core.id
# Add new thole nonbonded interaction for D-D, D-C, C-C for all existing
# drude types if this type is seen for the first time
if not type_drude in drude_dict:
# Bookkeeping of q, alphas and damping parameter
drude_dict[type_drude] = {}
drude_dict[type_drude]["q"] = q_drude
drude_dict[type_drude]["qc"] = p_core.q
drude_dict[type_drude]["alpha"] = alpha
drude_dict[type_drude]["thole_damping"] = thole_damping
drude_dict[type_drude]["core_type"] = p_core.type
# Save same information to get access to the parameters via core types
drude_dict[p_core.type] = {}
drude_dict[p_core.type]["q"] = -q_drude
drude_dict[p_core.type]["qc"] = p_core.q
drude_dict[p_core.type]["alpha"] = alpha
drude_dict[p_core.type]["thole_damping"] = thole_damping
drude_dict[p_core.type]["drude_type"] = type_drude
# Collect unique drude types
if not type_drude in drude_type_list:
drude_type_list.append(type_drude)
# Collect unique core types
if not p_core.type in core_type_list:
core_type_list.append(p_core.type)
# Collect unique drude ids
if not id_drude in drude_id_list:
drude_id_list.append(id_drude)
[docs]def add_thole_pair_damping(system, t1, t2, verbose=False):
"""
Calculates mixed Thole factors depending on Thole damping and polarization.
Adds non-bonded Thole interactions to the system.
Attributes
----------
system : Instance of :attr:`espressomd.System`
t1 : :obj:`int`
Type 1
t2 : :obj:`int`
Type 2
verbose : :obj:`bool`
Turns on verbosity.
"""
qq = drude_dict[t1]["q"] * drude_dict[t2]["q"]
s = 0.5 * (drude_dict[t1]["thole_damping"] + drude_dict[t2]["thole_damping"]
) / (drude_dict[t1]["alpha"] * drude_dict[t2]["alpha"])**(1.0 / 6.0)
system.non_bonded_inter[t1, t2].thole.set_params(scaling_coeff=s, q1q2=qq)
if verbose:
print("Added THOLE non-bonded interaction for types",
t1, "<->", t2, "S", s, "q1q2", qq)
[docs]def add_all_thole(system, verbose=False):
"""
Calls add_thole_pair_damping() for all necessary combinations to create the interactions.
Attributes
----------
system : Instance of :attr:`espressomd.System`
verbose : :obj:`bool`
Turns on verbosity.
"""
# drude <-> drude
for i in range(len(drude_type_list)):
for j in range(i, len(drude_type_list)):
add_thole_pair_damping(
system, drude_type_list[i], drude_type_list[j], verbose)
# core <-> core
for i in range(len(core_type_list)):
for j in range(i, len(core_type_list)):
add_thole_pair_damping(
system, core_type_list[i], core_type_list[j], verbose)
# drude <-> core
for i in drude_type_list:
for j in core_type_list:
add_thole_pair_damping(system, i, j, verbose)
[docs]def setup_and_add_drude_exclusion_bonds(system, verbose=False):
"""
Creates electrostatic short-range exclusion bonds for global exclusion
between Drude particles and core charges and adds the bonds to the cores.
Has to be called once after all Drude particles have been created.
Attributes
----------
system : Instance of :attr:`espressomd.System`
verbose: :obj:`bool`
Turns on verbosity.
"""
# All drude types need...
for td in drude_type_list:
#...exclusions with core
qd = drude_dict[td]["q"] # Drude charge
qc = drude_dict[td]["qc"] # Core charge
subtr_p3m_sr_bond = espressomd.interactions.BondedCoulombP3MSRBond(
q1q2=-qd * qc)
system.bonded_inter.add(subtr_p3m_sr_bond)
drude_dict[td]["subtr_p3m_sr_bonds_drude-core"] = subtr_p3m_sr_bond
if verbose:
print("Added drude-core p3m SR exclusion bond ",
subtr_p3m_sr_bond, "for drude", qd, "<-> core", qc, "to system")
for drude_id in drude_id_list:
core_id = core_id_from_drude_id[drude_id]
pd = system.part[drude_id]
pc = system.part[core_id]
bond = drude_dict[pd.type]["subtr_p3m_sr_bonds_drude-core"]
pc.add_bond((bond, drude_id))
if verbose:
print("Added drude-core p3m SR bond", bond,
"between ids", drude_id, "and", core_id)
[docs]def setup_intramol_exclusion_bonds(system, mol_drude_types, mol_core_types,
mol_core_partial_charges, verbose=False):
"""
Creates electrostatic short-range exclusion bonds for intramolecular exclusion
between Drude particles and partial charges of the cores. Has to be called once
after all Drude particles have been created.
Attributes
----------
system : Instance of :attr:`espressomd.System`
mol_drude_types : List of types of drude particles within the molecule
mol_core_types : List of types of core particles within the molecule
mol_core_partial_charges : List of partial charges of core particles within the molecule
verbose : :obj:`bool`
Turns on verbosity.
"""
# All drude types need...
for td in mol_drude_types:
drude_dict[td]["subtr_p3m_sr_bonds_intramol"] = {}
#...p3m sr exclusion bond with other partial core charges...
for tc, qp in zip(mol_core_types, mol_core_partial_charges):
#...excluding the drude core partner
if drude_dict[td]["core_type"] != tc:
qd = drude_dict[td]["q"] # Drude charge
subtr_p3m_sr_bond = espressomd.interactions.BondedCoulombP3MSRBond(
q1q2=-qd * qp)
system.bonded_inter.add(subtr_p3m_sr_bond)
drude_dict[td]["subtr_p3m_sr_bonds_intramol"][
tc] = subtr_p3m_sr_bond
if verbose:
print("Added intramolecular exclusion", subtr_p3m_sr_bond,
"for drude", qd, "<-> core", qp, "to system")
[docs]def add_intramol_exclusion_bonds(system, drude_ids, core_ids, verbose=False):
"""
Applies electrostatic short-range exclusion bonds for the given ids.
Has to be applied for all molecules.
Attributes
----------
system : Instance of :attr:`espressomd.System`
drude_ids : IDs of Drude particles within a molecule.
core_ids : IDs of core particles within a molecule.
verbose : :obj:`bool`
Turns on verbosity.
"""
for drude_id in drude_ids:
for core_id in core_ids:
if core_id_from_drude_id[drude_id] != core_id:
pd = system.part[drude_id]
pc = system.part[core_id]
bond = drude_dict[pd.type][
"subtr_p3m_sr_bonds_intramol"][pc.type]
pd.add_bond((bond, core_id))
if verbose:
print("Added subtr_p3m_sr bond", bond,
"between ids", drude_id, "and", core_id)