Saving average charge frequently in titration simulation

[udayardahal]

Dear All,

I am running a titration simulation using MARTINI-type model. I have a couple of questions:

1. How to implement explicitly the interaction if I don't want to use lorentz_berthelot mixing rule? Is there a way where I can define interaction between different beads like how we define for single atoms in atom list?

2. While running the titration simulation; I would like to see the instantaneous fluctuation of the charge groups and therefore I want to save a trajectory of average charges rather than a single output file at the end of the simulation. If it is possible, how do i do that?

I appreciate your answer and suggestion.

Thank you.

[mlund]

1. If I understand it correctly you want to specify the epsilon/sigma combination for a specific set of atom types? This is possible using the custom keyword, see for example this section in the manual.

2. Have a look at the chargefluctuations analysis, described here. This will give you detailed per-atom analysis and you can use the pqrfile keyword to save a PQR with average charges. To visualise charge fluctuations, you could also have a look at this analysis.

Do let us know if this worked for you, or if you need something else.

[udayardahal]

Thank you. I will update how it goes once I implement as you suggested .

[udayardahal]

Dear @mlund
I am able to get the charge fluctuations and implement the custom non-bonded interaction. Thank you for the help
I have another question:
When simulating a poly-electrolyte titration, I don't see much movement of the polyelectrolyte after few initial movements. After few strong movements it seems like the polymer has some sort of position restraint. Also, even when all beads of the polyelectrolyte is charged, I see a collapsed conformation of the particle. This may not be a "faunus" issue but I am trying to understand what could have been wrong in the simulation.
Also, if I want to use the implicit water, how do I do it? Does mentioning of epsr=80 is sufficient when we don't use explicit water?

My poly-electrolyte has 1-side chain beads for each backbone bead which is charged or neutral.
Here is my script for the simulation:
#!/usr/bin/env yason.py

This sets up a small, rigid molecule with titratable sites.

Titration is performed with implicit protons, coupled to

a grand canonical bath of 1:1 salt

{% set aNaCl = 0.15 %}
{% set pH = 1.0 %}

temperature: 300 # in Kelvin
random: {seed: hardware} # seed for random number generator
mcloop: {macro: 10, micro: 6000} # number of steps in MC simulation
geometry: {type: cuboid, length: [70.,70.,70.]}
atomlist:
# properties for all interaction sites, here single amino acids
- C1: {q: 0, sigma: 4.7, eps: 3.5, mw: 72, dp: 2}
- Qd: {q: 1, sigma: 4.7, eps: 5.0, mw: 72, dp: 2}
- Qd0: {q: 0, sigma: 4.7, eps: 5.0, mw: 72, dp: 2}
#- OW: {q: -0.8476, sigma: 3.166, eps: 0.650, mw: 15.999}
#- HW: {q: 0.4238, sigma: 2, eps: 0, mw: 1.007}
- OW: {q: 0, sigma: 4.7, eps: 4.0, mw: 24.0}
- HW1: {q: 0.46, sigma: 0, eps: 0, mw: 24.0}
- HW2: {q: -0.46, sigma: 0, eps: 0, mw: 24.0}
- na: {q: 1.0, sigma: 4.7, mw: 1, dp: 110, eps: 5.0}
- cl: {q: -1.0, sigma: 4.7, mw: 1, dp: 110, eps: 3.5}
- H+: { implicit: True, activity: {{ 10**(-pH) }} }

moleculelist:
- PAH:
structure:
- C1: [ 175.600, 167.900, 139.500]
- Qd: [ 175.000, 167.100, 137.100]
- C1: [ 173.800, 168.000, 141.600]
- Qd: [ 172.600, 166.000, 142.800]
- C1: [ 174.500, 171.000, 142.800]
- Qd: [ 173.600, 172.900, 141.300]
- C1: [ 177.100, 170.700, 144.600]
- Qd: [ 178.600, 168.600, 144.200]
- C1: [ 177.200, 173.000, 146.500]
- Qd: [ 175.700, 173.300, 148.600]
- C1: [ 180.300, 174.000, 147.400]
- Qd: [ 182.300, 172.600, 146.500]
- C1: [ 180.500, 175.900, 150.000]
- Qd: [ 178.800, 177.800, 150.600]
- C1: [ 183.000, 177.100, 151.500]
- Qd: [ 185.300, 176.600, 150.600]
- C1: [ 182.700, 176.800, 154.400]
- Qd: [ 184.300, 175.100, 155.500]
- C1: [ 181.300, 178.700, 156.100]
- Qd: [ 180.700, 178.800, 158.600]
- C1: [ 181.700, 181.900, 154.600]
- Qd: [ 184.000, 181.900, 153.300]
- C1: [ 179.000, 183.200, 155.200]
- Qd: [ 178.800, 184.000, 157.600]
- C1: [ 177.700, 185.100, 152.900]
- Qd: [ 179.500, 185.000, 151.100]
- C1: [ 175.500, 187.000, 152.800]
- Qd: [ 173.600, 187.100, 154.600]
- C1: [ 176.500, 190.200, 152.200]
- Qd: [ 175.900, 192.400, 153.400]
- C1: [ 178.700, 192.000, 150.000]
- Qd: [ 179.700, 194.100, 151.100]
- C1: [ 180.300, 192.000, 146.800]
- Qd: [ 182.800, 192.000, 147.700]
- C1: [ 181.300, 192.600, 143.200]
- Qd: [ 181.500, 195.200, 142.800]
- C1: [ 183.900, 190.700, 142.300]
- Qd: [ 183.900, 188.900, 144.200]
- C1: [ 187.000, 191.100, 141.100]
- Qd: [ 188.500, 193.200, 140.900]
bondlist:
- harmonic: {index: [ 0, 2], k: 0.27, req: 2.700}
- harmonic: {index: [ 2, 4], k: 0.27, req: 2.700}
- harmonic: {index: [ 4, 6], k: 0.27, req: 2.700}
- harmonic: {index: [ 6, 8], k: 0.27, req: 2.700}
- harmonic: {index: [ 8,10], k: 0.27, req: 2.700}
- harmonic: {index: [10,12], k: 0.27, req: 2.700}
- harmonic: {index: [12,14], k: 0.27, req: 2.700}
- harmonic: {index: [14,16], k: 0.27, req: 2.700}
- harmonic: {index: [16,18], k: 0.27, req: 2.700}
- harmonic: {index: [18,20], k: 0.27, req: 2.700}
- harmonic: {index: [20,22], k: 0.27, req: 2.700}
- harmonic: {index: [22,24], k: 0.27, req: 2.700}
- harmonic: {index: [24,26], k: 0.27, req: 2.700}
- harmonic: {index: [26,28], k: 0.27, req: 2.700}
- harmonic: {index: [28,30], k: 0.27, req: 2.700}
- harmonic: {index: [30,32], k: 0.27, req: 2.700}
- harmonic: {index: [32,34], k: 0.27, req: 2.700}
- harmonic: {index: [34,36], k: 0.27, req: 2.700}
- harmonic: {index: [36,38], k: 0.27, req: 2.700}
- harmonic: {index: [ 0, 1], k: 100.0, req: 2.700}
- harmonic: {index: [ 2, 3], k: 100.0, req: 2.700}
- harmonic: {index: [ 4, 5], k: 100.0, req: 2.700}
- harmonic: {index: [ 6, 7], k: 100.0, req: 2.700}
- harmonic: {index: [ 8, 9], k: 100.0, req: 2.700}
- harmonic: {index: [10,11], k: 100.0, req: 2.700}
- harmonic: {index: [12,13], k: 100.0, req: 2.700}
- harmonic: {index: [14,15], k: 100.0, req: 2.700}
- harmonic: {index: [16,17], k: 100.0, req: 2.700}
- harmonic: {index: [18,19], k: 100.0, req: 2.700}
- harmonic: {index: [20,21], k: 100.0, req: 2.700}
- harmonic: {index: [22,23], k: 100.0, req: 2.700}
- harmonic: {index: [24,25], k: 100.0, req: 2.700}
- harmonic: {index: [26,27], k: 100.0, req: 2.700}
- harmonic: {index: [28,29], k: 100.0, req: 2.700}
- harmonic: {index: [30,31], k: 100.0, req: 2.700}
- harmonic: {index: [32,33], k: 100.0, req: 2.700}
- harmonic: {index: [34,35], k: 100.0, req: 2.700}
- harmonic: {index: [36,37], k: 100.0, req: 2.700}
- harmonic: {index: [38,39], k: 100.0, req: 2.700}
excluded_neighbours: 1
- water:
structure:
- OW: [1.060, 1.015, 0.979]
- HW1: [0.964, 0.921, 0.941]
- HW2: [0.975, 1.063, 1.080]
# - OW: [2.30, 6.28, 1.13]
# - HW1: [ 1.37,6.26, 1.50]
# - HW2: [2.31, 5.89, 0.21]
bondlist:
- harmonic: {index: [0,1], k: 100, req: 0.14}
- harmonic: {index: [0,2], k: 100, req: 0.14}
- Na+: {atoms: [na], atomic: true, activity: {{ aNaCl }} }
- Cl-: {atoms: [cl], atomic: true, activity: {{ aNaCl }} }

insertmolecules:
- PAH: {N: 1}
- water: {N: 100}
- Na+: {N: 30}
- Cl-: {N: 50}

reactionlist:
# H+ is implicit and to maintain electroneutrality
# we accompany swap moves with insertion/deletion
# of either a chloride or a sodium ion.
# Effectively this means that we can interconvert
# H+ <-> Na+ <-> Cl-. Excess salt pairs are removed
# by the final NaCl GC reaction.
#
# pKa values from dx.doi.org/10.1110/ps.051840806
#
#- Qd = Na+ + Qd0 + H+: {pK: 10.23}
#- Qd + Cl- = Qd0 + H+: {pK: 10.23}
- Qd = Qd0 + H+: {pK: 10.23} # not electroneutral
- = Na+ + Cl-: {} # activates grand canonical salt

energy:
- nonbonded:
default:
- coulomb: {epsr: 1.3, type: ewald, cutoff: 14, debyelength: 100, alpha: 0.2, kcutoff: 7 }
#- coulomb: {epsr: 1, type: plain}
- lennardjones:
mixing: LB
custom:
- OW Qd: {eps: 5.0, sigma: 4.7}
- OW Qd0: {eps: 5.0, sigma: 4.7}
- OW na: {eps: 5.0, sigma: 4.7}
- OW cl: {eps: 5.0, sigma: 4.7}
- na cl: {eps: 4.0, sigma: 4.7}
- Qd cl: {eps: 4.0, sigma: 4.7}
- Qd0 cl: {eps: 4.0, sigma: 4.7}
- C1 Qd: {eps: 2.5, sigma: 4.7}
- C1 Qd0: {eps: 2.5, sigma: 4.7}
- C1 na: {eps: 2.5, sigma: 4.7}
- C1 cl: {eps: 2.5, sigma: 4.7}
- isobaric: {P/atm: 1}
- bonded: {}

moves:
- rcmc: {repeat: 10} # reactive ensemble for proton titration
- crankshaft: {molecule: PAH, dprot: 0.5, repeat: 1}
- moltransrot: {molecule: PAH, dp: 20, dprot: 0.5, repeat: 1}
- pivot: {molecule: PAH, dprot: 0.5, repeat: 1}
- transrot: {molecule: PAH, repeat: 1}
- transrot: {molecule: Na+, repeat: N} # move individual salt beads
- transrot: {molecule: Cl-, repeat: N} # move individual salt beads
- moltransrot: {molecule: water, dp: 2, dprot: 2, repeat: N}

analysis:
- savestate: {file: confout.pqr}
- sanity: {nstep: 10} # is the simulation OK?
- chargefluctuations: {nstep: 10, molecule: PAH, pqrfile: avgcharge.pqr} # charge analysis of all sites
- multipole: {nstep: 5} # molecular charge, capacitance and dipole moment
- systemenergy: {file: energy_3.dat, nstep: 100} # track total system energy
- density: {nstep: 5, nskip: 10} # monitors salt density
- savestate: {file: state_3.json}
- xtcfile: {file: traj_3.xtc, nstep: 10}
- qrfile: {file: qrtraj_3.dat, nstep: 10}

[mlund]

The simulation setup seems to be with explicit solvent. Is this what you want, or do you want implicit solvent?

• If implicit, then

  1. remove water from insert_molecules
  2. use epsr: 80 in nonbonded. You could also consider using plain instead of ewald which will be faster. Generally, Ewald is not needed in the implicit setup you have.
  3. Your LJ epsilon values of ~5 kJ/mol seem large for implicit solvent and would likely cause a collapsed chain. As a rule of thumb vdW interactions are ~10 times lower in implicit solvent than in explicit solvent (if we're talking about water).
  4. remove the isobaric energy

• if explicit solvent the setup is problematic in that the MC trial energies for both titration and GC salt insertion/deletion are on the order of Born energies (many tens of kT) and therefore almost impossible to sample.

[mlund]

Please note #282 which could be useful for simulations of MARTINI models combined with MC titration. @himalayanpride can I close this issue?

[udayardahal]

Please note #282 which could be useful for simulations of MARTINI models combined with MC titration. @himalayanpride can I close this issue?

Sorry for late reply. You can close this issue.