This module is for performing simulations (e.g. as part of Thermodynamic Integration) where some interactions are a linear function of a parameter \(\lambda\).

\[U(\lambda) = (1-\lambda)U_C^A\]

where \(U_C^A\) is the standard Reaction Field interaction. This allows one to perform TI where the charges in TI state A (\(\lambda=0\)) are the particle charges contained in the particle property charge and the charges in TI state B (\(\lambda=1\)) are zero.

The user specifies a list of particles, pidlist. For all pairs of particles with particletypes interacting via this potential, the RF interaction between two particles i and j is calculated as follows:

if (i not in pidlist) and (j not in pidlist):
\(U_{RF}\) (full RF interaction)
if (i in pidlist) and (j in pidlist):
if annihilate==True:
\((1-\lambda)U_{RF}\) (RF interaction scaled by 1-lambda)
if annihilate==False:
\(U_{RF}\) (full RF interaction)
if (i in pidlist) xor (j in pidlist):
\((1-\lambda)U_{RF}\) (RF interaction scaled by 1-lambda)

The default is annihilation (completely turning off charges of particles in pidlist in state B, so that interactions within pidlist are turned off and also cross-interactions between particles in pidlist and particles in the rest of the system). The alternative is decoupling (only cross-interactions between particles in pidlist and particles in the rest of the system are turned off. Interactions within pidlist are not affected.) If annihilation==False, then decoupling is performed. See:

Exclusions apply as normal, i.e. interactions are only calculated for pairs of particles not already excluded.

So far only VerletListAdressReactionFieldGeneralizedTI is implemented, however VerletListReactionFieldGeneralizedTI, VerletListHadressReactionFieldGeneralizedTI, etc. can also be easily implemented.

The \(\lambda\) (lambdaTI) parameter used here should not be confused with the \(\lambda\) (lambda_adr) particle property used in AdResS simulations.

See also the Thermodynamic Integration tutorial.

Example python script:

>>> #value of lambda
>>> lambdaTI = 0.3
>>> #construct RF potential with parameters prefactor,kappa,epsilon1,epsilon2,cutoff as in standard RF interaction
>>> pot = espressopp.interaction.ReactionFieldGeneralizedTI(prefactor=prefactor, kappa=kappa, epsilon1=epsilon1, epsilon2=epsilon2, cutoff=rc, lambdaTI=lambdaTI, annihilate=False)
>>> #add list of indices of particles whose charge is 0 in TI state B
>>> pidlist = [1,2,3,4]
>>> pot.addPids(pidlist)
>>> #create interaction using VerletListAdress object and FixedTupleListAdress object
>>> qq_adres_interaction=espressopp.interaction.VerletListAdressReactionFieldGeneralizedTI(verletlist, ftpl)
>>> #loop over list of all types for particles interacting with this atomistic potential
>>> for i in types:
>>>   for k in types:
>>>     qq_adres_interaction.setPotentialAT(type1=i, type2=k, potential=pot)
>>> system.addInteraction(qq_adres_interaction)

During the MD run, one can then calculate the derivative of the RF energy wrt lambda

>>> #calculate dU/dlambda
>>> dUdl = qq_adres_interaction.computeEnergyDeriv()
espressopppp.interaction.ReactionFieldGeneralizedTI(prefactor, kappa, epsilon1, epsilon2, cutoff, lambdaTI, annihilate)
  • prefactor (real) – (default: 1.0) RF parameter
  • kappa (real) – (default: 0.0) RF parameter
  • epsilon1 (real) – (default: 1.0) RF parameter
  • epsilon2 (real) – (default: 80.0) RF parameter
  • cutoff (real) – (default: infinity) interaction cutoff
  • lambdaTI (real) – (default: 0.0) TI lambda parameter
  • annihilate (bool) – (default: True) switch between annihilation and decoupling
Parameters:pidlist (python list) – list of particle ids of particles whose charge is zero in state B
espressopppp.interaction.VerletListAdressReactionFieldGeneralized(vl, fixedtupleList)
  • vl (VerletListAdress object) – Verlet list
  • fixedtupleList (FixedTupleListAdress object) – list of tuples describing mapping between CG and AT particles
espressopppp.interaction.VerletListAdressReactionFieldGeneralized.setPotentialAT(type1, type2, potential)
  • type1 (int) – atomtype
  • type2 (int) – atomtype
  • potential (Potential) – espressopppp potential
espressopppp.interaction.VerletListAdressReactionFieldGeneralized.setPotentialCG(type1, type2, potential)
  • type1 (int) – atomtype
  • type2 (int) – atomtype
  • potential (Potential) – espressopppp potential
class espressopp.interaction.ReactionFieldGeneralizedTI.ReactionFieldGeneralizedTI

The ReactionFieldGeneralizedTI potential.