[gmx-users] calculation of pKa
Rui Neves
rppneves at gmail.com
Mon Aug 29 16:31:18 CEST 2016
Dear Gromacs users,
I am trying to perform a pKa calculation of a Cys in a thioredoxin-protein
(Trx) using the FEP implementation of Gromacs 5.1.0.
I am transforming the charge in CYS to that of a CYM in 10 equally-spaced
lambda steps, and then turning off the vdw radii of the hydrogen (HS) in
the thiol terminus of Cys in 20 equally-spaced lambda steps.
I have had no problems concerning the anhilation and redistribution of
charge in Cys.
However, I have had some problems with the anhilation of the vdw radii of
HS.
In this process I turn a HS into a dummy (DU) with no LJ energy of
interaction.
During the BAR analysis with the gmx bar program, I frequently come across
this warning:
*WARNING: Some of these results violate the Second Law of Thermodynamics: *
* This is can be the result of severe undersampling, or (more
likely)*
* there is something wrong with the simulations.*
I have tried to read about this in forums. I often read that it is related
to a negative entropy for states A and B to 'interact' with one another,
but there is no simple answer to solve this problem.
I have revised my mdp settings and have managed to solve the issue for the
case of the system comprising only Cys; however, I still cannot solve this
issue for the Trx system.
I tried to extend the simulation from 1 ns to 4 ns/lambda, but I keep
getting the same warning (I attach the bar plot to this mail).
*Temperature: 300 K*
*Detailed results in kT (see help for explanation):*
* lam_A lam_B DG +/- s_A +/- s_B +/-
stdev +/- *
* 10 11 -0.0359 0.0005 0.0060 0.0018 0.0070 0.0018
0.1060 0.0006*
* 11 12 -0.0427 0.0004 0.0120 0.0010 0.0130 0.0010
0.1103 0.0007*
* 12 13 -0.0414 0.0007 0.0003 0.0005 0.0012 0.0006
0.1032 0.0007*
* 13 14 -0.0384 0.0008 0.0080 0.0005 0.0087 0.0005
0.0986 0.0008*
* 14 15 -0.0394 0.0004 0.0053 0.0005 0.0059 0.0005
0.0959 0.0004*
* 15 16 -0.0322 0.0009 -0.0004 0.0007 -0.0001 0.0007
0.0834 0.0007*
* 16 17 -0.0267 0.0004 0.0046 0.0009 0.0047 0.0009
0.0708 0.0002*
* 17 18 -0.0214 0.0011 -0.0002 0.0012 -0.0001 0.0012
0.0598 0.0006*
* 18 19 -0.0162 0.0003 0.0025 0.0016 0.0026 0.0016
0.0510 0.0006*
* 19 20 -0.0124 0.0006 0.0007 0.0010 0.0008 0.0010
0.0424 0.0004*
* 20 21 -0.0094 0.0006 0.0022 0.0007 0.0022 0.0007
0.0360 0.0004*
* 21 22 -0.0046 0.0007 -0.0014 0.0003 -0.0014 0.0003
0.0300 0.0004*
* 22 23 0.0010 0.0005 -0.0002 0.0007 -0.0002 0.0007
0.0239 0.0002*
* 23 25 0.0083 0.0012 0.0031 0.0007 0.0019 0.0004
0.0364 0.0002*
* 25 26 0.0071 0.0003 0.0002 0.0002 0.0003 0.0004
0.0138 0.0001*
* 26 27 0.0093 0.0002 -0.0004 0.0002 -0.0004 0.0002
0.0109 0.0001*
* 27 28 0.0115 0.0001 -0.0002 0.0002 -0.0002 0.0002
0.0090 0.0001*
* 28 29 0.0126 0.0001 0.0004 0.0002 0.0004 0.0001
0.0077 0.0001*
* 29 30 0.0130 0.0001 0.0004 0.0000 0.0004 0.0000
0.0065 0.0000*
*WARNING: Some of these results violate the Second Law of Thermodynamics: *
* This is can be the result of severe undersampling, or (more
likely)*
* there is something wrong with the simulations.*
*Final results in kJ/mol:*
*point 10 - 11, DG -0.0895 +/- 0.0013*
*point 11 - 12, DG -0.1066 +/- 0.0010*
*point 12 - 13, DG -0.1032 +/- 0.0019*
*point 13 - 14, DG -0.0959 +/- 0.0019*
*point 14 - 15, DG -0.0983 +/- 0.0009*
*point 15 - 16, DG -0.0802 +/- 0.0022*
*point 16 - 17, DG -0.0666 +/- 0.0010*
*point 17 - 18, DG -0.0533 +/- 0.0028*
*point 18 - 19, DG -0.0403 +/- 0.0007*
*point 19 - 20, DG -0.0309 +/- 0.0014*
*point 20 - 21, DG -0.0235 +/- 0.0016*
*point 21 - 22, DG -0.0115 +/- 0.0017*
*point 22 - 23, DG 0.0025 +/- 0.0013*
*point 23 - 25, DG 0.0206 +/- 0.0031*
*point 25 - 26, DG 0.0176 +/- 0.0008*
*point 26 - 27, DG 0.0232 +/- 0.0004*
*point 27 - 28, DG 0.0288 +/- 0.0002*
*point 28 - 29, DG 0.0314 +/- 0.0002*
*point 29 - 30, DG 0.0325 +/- 0.0001*
*total 10 - 30, DG -0.6432 +/- 0.0037*
I have also analysed different timeblocks of the simulation, to no avail.
I have also tested the inclusion of the couple-moltype keyword, but the
problem remains. I get different free-energy estimates though. Does anyone
know what exactly means the couple-moltype keyword?
My mdp settings are as follow:
*; Run control*
*integrator = sd ; Langevin dynamics*
*tinit = 0*
*dt = 0.002*
*nsteps = 2000000 ; 4 ns*
*nstcomm = 100*
*comm-mode = Linear*
*; Output control*
*nstxout = 0 *
*nstvout = 0 *
*nstfout = 0*
*nstlog = 1000*
*nstenergy = 1000*
*nstxout-compressed = 1000*
*nstcalcenergy = 100*
*; Neighborsearching and short-range nonbonded interactions*
*cutoff-scheme = verlet*
*nstlist = 10*
*ns_type = grid*
*pbc = xyz*
*rlist = 1.0*
*; Electrostatics*
*coulombtype = PME*
*rcoulomb-switch = 0.96*
*rcoulomb = 1.0*
*coulomb-modifier = Potential-shift-Verlet*
*; van der Waals*
*vdwtype = cutoff*
*vdw-modifier = Potential-shift-Verlet*
*rvdw-switch = 0.96*
*rvdw = 1.0*
*; Apply long range dispersion corrections for Energy and Pressure*
*DispCorr = EnerPres*
*; Spacing for the PME/PPPM FFT grid*
*fourierspacing = 0.12*
*; EWALD/PME/PPPM parameters*
*pme_order = 6*
*ewald_rtol = 1e-06*
*epsilon_surface = 0*
*; Temperature coupling*
*; tcoupl is implicitly handled by the sd integrator*
*tc_grps = Protein Water_and_ions*
*tau_t = 1.0 1.0*
*ref_t = 300 300*
*; Pressure coupling is on for NPT*
*Pcoupl = Parrinello-Rahman *
*tau_p = 1.0*
*compressibility = 4.5e-05*
*ref_p = 1.0 *
*; Free energy control stuff*
*free_energy = yes*
*init_lambda_state = $LAMBDA*
*delta_lambda = 0*
*calc_lambda_neighbors = -1 ; all neighboring windows*
*; Vectors of lambda specified here*
*; Each combination is an index that is retrieved from init_lambda_state
for each simulation*
*; init_lambda_state 0 1 2 3 4 5 6 7 8 9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
25 26 27 28 29 30*
*vdw_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65
0.70 0.75 0.80 0.85 0.90 0.95 1.00*
*coul_lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80
0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00*
*; We are not transforming any bonded or restrained interactions*
*bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00*
*restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00*
*; Masses are not changing (particle identities are the same at lambda = 0
and lambda = 1)*
*mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00*
*; Not doing simulated temperting here*
*temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00*
*; Options for the decoupling*
*sc-alpha = 0.5*
*sc-coul = no ; linear interpolation of Coulomb
(none in this case)*
*sc-power = 1.0*
*sc-sigma = 0.3 *
*couple-moltype = Protein_chain_A*
*couple-intramol = yes*
*nstdhdl = 10*
*; Do not generate velocities*
*gen_vel = no *
*; options for bonds*
*constraints = h-bonds *
*; Type of constraint algorithm*
*constraint-algorithm = lincs*
*; Constrain the starting configuration*
*; since we are continuing from NPT*
*continuation = yes *
*; Highest order in the expansion of the constraint coupling matrix*
*lincs-order = 12*
I also post the changes that I made to create the B-state for the Cys32
residue in Trx.
*; Include forcefield parameters*
*#include "amber99sb-ildn.ff/forcefield.itp"*
*[ atomtypes ]*
*;name bond_type mass charge ptype sigma epsilon
Amb*
* DU HS 0.00000 0.00000 A 1.06908e-01 0.00000e+00 ;
0.60 0.0000*
*[ moleculetype ]*
*; Name nrexcl*
*Protein_chain_A 3*
*[ atoms ]*
* ( ... )*
*; residue 32 CYS rtp CYS q 0.0*
* 478 N 32 CYS N 478 -0.4157 14.01 N
-0.4157 14.01 ; qtot -1.416*
* 479 H 32 CYS H 479 0.2719 1.008 H
0.2719 1.008 ; qtot -1.144*
* 480 CT 32 CYS CA 480 0.0213 12.01 CT
-0.0351 12.01 ; qtot -1.123*
* 481 H1 32 CYS HA 481 0.1124 1.008 H1
0.0508 1.008 ; qtot -1.01*
* 482 CT 32 CYS CB 482 -0.1231 12.01 CT
-0.2413 12.01 ; qtot -1.133*
* 483 H1 32 CYS HB1 483 0.1112 1.008 H1
0.1122 1.008 ; qtot -1.022*
* 484 H1 32 CYS HB2 484 0.1112 1.008 H1
0.1122 1.008 ; qtot -0.9108*
* 485 SH 32 CYS SG 485 -0.3119 32.06 SH
-0.8844 32.06 ; qtot -1.223*
* 486 HS 32 CYS HG 486 0.1933 1.008 DU
0.0000 1.008 ; qtot -1.029*
* 487 C 32 CYS C 487 0.5973 12.01 C
0.5973 12.01 ; qtot -0.4321*
* 488 O 32 CYS O 488 -0.5679 16 O
-0.5679 16 ; qtot -1*
* (...)*
Do you have any ideas on what else I can try?
Thank you so much for your time.
Regards,
Rui PP Neves
More information about the gromacs.org_gmx-users
mailing list