Harmonic bond breaking on Serine phosphorylation reaction

[directysj]

Hello Faunus team,

In our system, we would like to phosphorylate Serine, therefore, first I tested it in simplistic system to check how it could behave. I am giving the tested possibilities below:

1). Single serine with counter ions. 2). Multiple non-bonded serines with counter ions. 3). Two serines covalently bonded as a dumbbell for simplicity.

Unfortunately, the bond is breaking in the third case. I have tested many possibilities by varying average bond length, diameter, and molecular weight, switching on/off non-bonded interactions, to see which of them is creating the problem. Surprisingly, it pointed out that when molecular weights of SER and SEP are different (while covalently bonded), the bond is broken and the simulation is ending by printing an error message. I have looked the configuration of the system in VMD which shows that the bond stretches too much. If molecular weights of SER and SEP are equal then the simulation is running nicely. I am attaching the JSON script below for your reference.

#!/usr/bin/env yason.py
{% set eps = 0.10 %} # interaction energy in kJ/mol
{% set epsSalt = 0.10 %} # interaction energy in kJ/mol
{% set pH = 7.0 %}
{% set pOH = 14 - pH %}
{% set phosActivity = 0.01 %} # activity in mol/l
{% set h2oActivity = 0.0000001 %} # activity in mol/l
{% set saltActivity = 0.154 %} # activity in mol/l
{% set debyelength = 3.04/saltActivity**0.5 %}
{% set dpsalt = 20 %} # displacement in angstrom 
{% set dpSer = 1.0 %}
{% set Lx = 100.0 %}
{% set Lz = 100.0 %}

temperature: 300 
random: {seed: hardware}
geometry: {type: cuboid, length: [{{ Lx }}, {{ Lx }}, {{ Lz }}]}
mcloop: {macro: 100, micro: 1000}

atomlist:
    - NTR:     { q: 1,  sigma: 4.0, mw: 14,  dp: {{ dpSer }},  eps: {{ eps }} }
    - CTR:     { q: -1, sigma: 4.0, mw: 16,  dp: {{ dpSer }},  eps: {{ eps }} }
    - SER:     { q: 0,  sigma: 6.6, mw: 82,  dp: {{ dpSer }},  eps: {{ eps }} }
    - TYR:     { q: -1, sigma: 8.2, mw: 154, dp: {{ dpSer }},  eps: {{ eps }} }
    - SEP:     { q: -2, sigma: 8.4, mw: 161, dp: {{ dpSer }},  eps: {{ eps }} }
    - HSEP:    { q: -1, sigma: 8.4, mw: 161, dp: {{ dpSer }},  eps: {{ eps }} }
    - H2SEP:   { q: 0,  sigma: 8.4, mw: 161, dp: {{ dpSer }},  eps: {{ eps }} }
    - na:      { q: 1,  sigma: 4.6, mw: 11,  dp: {{ dpsalt }}, eps: {{ epsSalt }} }
    - cl:      { q: -1, sigma: 4.6, mw: 35,  dp: {{ dpsalt }}, eps: {{ epsSalt }} }
    - H3PO4:   { q: 0,  implicit: True, activity: {{ phosActivity }} }
    - H2PO4:   { q: -1, implicit: True, activity: {{ phosActivity }} }
    - HPO4:    { q: -2, implicit: True, activity: {{ phosActivity }} }
    - H2O:     { q: 0,  implicit: True, activity: {{ h2oActivity }} }
    - H+:      { q: 1,  implicit: True, activity: {{ 10**(-pH) }} }
    - OH-:     { q: -1, implicit: True, activity: {{ 10**(-pOH) }} }

moleculelist:
    - NaCl:  {atoms: [na, cl], atomic: True, activity: {{ saltActivity }} }
    - Serine:
         structure: {fasta: "nSSc", k: 10.0, req: 8.0}

reactionlist:
    # H+ is implicit and used for the titration moves
    - "H3PO4 = H2PO4 + H+":      {pK: 2.12}
    - "H2PO4 = HPO4 + H+":       {pK: 7.21}
    - "H2O = H+ + OH-":          {pK: 14.0}
    # Phosphorylation of Serine 
    - "SER + HPO4 = SEP + H2O":  {pK: 5.78} 
    - "HSEP = SEP + H+":         {pK: 5.78}
    - "H2SEP = HSEP + H+":       {pK: 2.19}

insertmolecules:
    - Serine: {N: 1} # Serine has zero charge
    - NaCl: {N: 100} # Salt

energy:
    - bonded: {}
    - nonbonded_coulomblj:
                openmp: [i2all]
                lennardjones: {mixing: LB}
                coulomb: {epsr: 78.7, type: yukawa, shift: true, debyelength: {{ debyelength }}, cutoff: {{ 3.0*debyelength }} }
moves:
    - rcmc:         {}
    - transrot:     {molecule: Serine, repeat: N}
    - transrot:     {molecule: NaCl, repeat: N}

analysis:
    - sanity: {nstep: 1}
    - multipole: {nstep: 500}
    - reactioncoordinate: {file: systemcharge.dat, type: system, property: Q, nstep: 500}
    - xtcfile: {file: traj.xtc, nstep: 10000}
    - savestate: {file: confout.pqr}
    - savestate: {file: state.json, nstep: 50000}

Can I ask your attention/help in this regard?

Sincerely,

[mlund]

Hmm, not quite sure where it fails. A few things:

1. Could you reduce the input to a bare minimum where you still see the issue?
2. One possible source of error could be that the mass center of Serine is not updated during the swap move. Do you still see it if SER and SEP have equal mw? Confirmed that it's the mass center.
3. Why include water auto-dissociation? I suspect water is in great excess. How was the water activity selected? I suggest removing reactions with H2O and OH- for now.

Update: removing the fourth reaction solves the issue, hinting that the mass center is not updated during this move. This needs to be fixed, but in your case, I don't see why you would need that reaction, see (3). That is, just removing "H2O" from the fourth reaction should fix the issue for you.

[directysj]

Thank you so much! Yes, it is the mass center problem. I am not getting error after deletion of water. Actually, I found this chemical reaction for the phosphorylation of SER and Faunus had no restrictions for reactants and products, that's why I kept in this way.

[mlund]

:-) Please see above PR that fixes the mass center problem. If OK, I'll merge it into the master branch.

[directysj]

Now, the problem is that the reaction is not taking place. See the section of out.json. But NO ERROR!!

"reactions": {
           "SER + HPO4 = SEP": {
             "acceptance -->": 0.0,
             "acceptance <--": null,
             "attempts": 50053
           }
}

[mlund]

Isn't that because particle diameter increases (SER -> SEP) so that a large overlap energy is associated with a forward reaction? The choice of pK and activity will also make the reaction virtually impossible. Testing with e.g. - "SEP = SER + HPO4": {pK: 5.7} and a HPO4 activity of 10**(-5.7) gives a fifty/fifty distribution as expected

[directysj]

To avoid overlapping, I am keeping initial interatomic distance between SER atoms in such a way that they should not overlap if the forward reaction takes place. Yes, I will have a look at the activities of the implicit reactants and products.