Source code for espressomd.drude_helpers

# Copyright (C) 2010-2022 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 . import interactions
from .code_features import has_features


[docs] class DrudeHelpers: """ Provides helper functions to assist in the creation of Drude particles. """ def __init__(self): # Dict with Drude type information self.drude_dict = {} # Lists with unique Drude and core types self.core_type_list = [] self.drude_type_list = [] # Get core id from Drude id self.core_id_from_drude_id = {} # Drude IDs self.drude_id_list = []
[docs] def add_drude_particle_to_core(self, system, harmonic_bond, thermalized_bond, p_core, type_drude, alpha, mass_drude, coulomb_prefactor, thole_damping=2.6, verbose=False): """ Adds a Drude particle with specified type and mass to the system and returns its `espressomd.particle_data.ParticleHandle`. Checks if different Drude particles have different types. Collects types/charges/polarizations/Thole factors for intramolecular core-Drude short-range exclusion and Thole interaction. Parameters ---------- system : :class:`espressomd.system.System` harmonic_bond: :class:`espressomd.interactions.HarmonicBond` Add this harmonic bond to between Drude particle and core thermalized_bond: :class:`espressomd.interactions.ThermalizedBond` Add this thermalized_bond to between Drude particle and core p_core: :class:`espressomd.particle_data.ParticleHandle` The existing core particle 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 :meth:`add_all_thole()` method is used. verbose : :obj:`bool` Turns on verbosity. Returns ------- :class:`espressomd.particle_data.ParticleHandle` The created Drude Particle. """ k = harmonic_bond.params["k"] q_drude = -1.0 * pow(k * alpha / coulomb_prefactor, 0.5) if has_features("THERMOSTAT_PER_PARTICLE"): if has_features("PARTICLE_ANISOTROPY"): gamma_off = [0.0, 0.0, 0.0] else: gamma_off = 0.0 drude_part = system.part.add(pos=p_core.pos, type=type_drude, q=q_drude, mass=mass_drude) if has_features("THERMOSTAT_PER_PARTICLE"): drude_part.gamma = gamma_off id_drude = drude_part.id if verbose: print( f"Adding to core {p_core.id} 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)) if has_features("THERMOSTAT_PER_PARTICLE"): p_core.gamma = gamma_off if type_drude in self.drude_dict and not ( self.drude_dict[type_drude]["q"] == q_drude and self.drude_dict[type_drude]["thole_damping"] == thole_damping): raise Exception( "Drude particles with different drude charges have to have different types for THOLE") self.core_id_from_drude_id[id_drude] = p_core.id # Add new Thole non-bonded interaction for D-D, D-C, C-C for all existing # Drude types if this type is seen for the first time if type_drude not in self.drude_dict: # Bookkeeping of q, alphas and damping parameter self.drude_dict[type_drude] = {} self.drude_dict[type_drude]["q"] = q_drude self.drude_dict[type_drude]["qc"] = p_core.q self.drude_dict[type_drude]["alpha"] = alpha self.drude_dict[type_drude]["thole_damping"] = thole_damping self.drude_dict[type_drude]["core_type"] = p_core.type # Save same information to get access to the parameters via core # types self.drude_dict[p_core.type] = {} self.drude_dict[p_core.type]["q"] = -q_drude self.drude_dict[p_core.type]["qc"] = p_core.q self.drude_dict[p_core.type]["alpha"] = alpha self.drude_dict[p_core.type]["thole_damping"] = thole_damping self.drude_dict[p_core.type]["drude_type"] = type_drude # Collect unique Drude types if type_drude not in self.drude_type_list: self.drude_type_list.append(type_drude) # Collect unique core types if p_core.type not in self.core_type_list: self.core_type_list.append(p_core.type) # Collect unique Drude ids if id_drude not in self.drude_id_list: self.drude_id_list.append(id_drude) return drude_part
[docs] def add_thole_pair_damping(self, system, t1, t2, verbose=False): """ Calculates mixed Thole factors depending on Thole damping and polarization. Adds non-bonded Thole interactions to the system. Parameters ---------- system : :class:`espressomd.system.System` t1 : :obj:`int` Type 1 t2 : :obj:`int` Type 2 verbose : :obj:`bool` Turns on verbosity. """ qq = self.drude_dict[t1]["q"] * self.drude_dict[t2]["q"] s = 0.5 * (self.drude_dict[t1]["thole_damping"] + self.drude_dict[t2]["thole_damping"] ) / (self.drude_dict[t1]["alpha"] * self.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( f"Added THOLE non-bonded interaction for types {t1} <-> {t2} S {s} q1q2 {qq}")
[docs] def add_all_thole(self, system, verbose=False): """ Calls :meth:`add_thole_pair_damping()` for all necessary combinations to create the interactions. Parameters ---------- system : :class:`espressomd.system.System` verbose : :obj:`bool` Turns on verbosity. """ # Drude <-> Drude for i in range(len(self.drude_type_list)): for j in range(i, len(self.drude_type_list)): self.add_thole_pair_damping( system, self.drude_type_list[i], self.drude_type_list[j], verbose) # core <-> core for i in range(len(self.core_type_list)): for j in range(i, len(self.core_type_list)): self.add_thole_pair_damping( system, self.core_type_list[i], self.core_type_list[j], verbose) # Drude <-> core for i in self.drude_type_list: for j in self.core_type_list: self.add_thole_pair_damping(system, i, j, verbose)
[docs] def setup_and_add_drude_exclusion_bonds(self, 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. Parameters ---------- system : :class:`espressomd.system.System` verbose: :obj:`bool` Turns on verbosity. """ # All Drude types need... for td in self.drude_type_list: # ...exclusions with core qd = self.drude_dict[td]["q"] # Drude charge qc = self.drude_dict[td]["qc"] # Core charge subtr_sr_bond = interactions.BondedCoulombSRBond( q1q2=-qd * qc) system.bonded_inter.add(subtr_sr_bond) self.drude_dict[td]["subtr_sr_bonds_drude-core"] = subtr_sr_bond if verbose: print( f"Added drude-core SR exclusion bond {subtr_sr_bond} for drude {qd} <-> core {qc} to system") for drude_id in self.drude_id_list: core_id = self.core_id_from_drude_id[drude_id] pd = system.part.by_id(drude_id) pc = system.part.by_id(core_id) bond = self.drude_dict[pd.type]["subtr_sr_bonds_drude-core"] pc.add_bond((bond, drude_id)) if verbose: print( f"Added drude-core SR bond {bond} between ids {drude_id} and {core_id}")
[docs] def setup_intramol_exclusion_bonds(self, 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. Parameters ---------- system : :class:`espressomd.system.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: self.drude_dict[td]["subtr_sr_bonds_intramol"] = {} # ... 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 self.drude_dict[td]["core_type"] != tc: qd = self.drude_dict[td]["q"] # Drude charge subtr_sr_bond = interactions.BondedCoulombSRBond( q1q2=-qd * qp) system.bonded_inter.add(subtr_sr_bond) self.drude_dict[td]["subtr_sr_bonds_intramol"][ tc] = subtr_sr_bond if verbose: print( f"Added intramolecular exclusion {subtr_sr_bond} for drude {qd} <-> core {qp} to system")
[docs] def add_intramol_exclusion_bonds( self, drude_parts, core_parts, verbose=False): """ Applies electrostatic short-range exclusion bonds for the given ids. Has to be applied for all molecules. Parameters ---------- system : :class:`espressomd.system.System` drude_parts : List of Drude particles within a molecule. core_parts : List of core particles within a molecule. verbose : :obj:`bool` Turns on verbosity. """ for drude_part in drude_parts: drude_id = drude_part.id for core_part in core_parts: core_id = core_part.id if self.core_id_from_drude_id[drude_id] != core_id: pd = drude_part pc = core_part bond = self.drude_dict[pd.type][ "subtr_sr_bonds_intramol"][pc.type] pd.add_bond((bond, core_id)) if verbose: print( f"Added subtr_sr bond {bond} between ids {drude_id} and {core_id}")