Solvation free energy calculation of organic molecule in binary mixtures and convergence

User discussions
1

Dear all,

I am trying to calculate the solvation free energy for an organic molecule (merocyanine) in pure water and water/acetonitrile binary mixture.

I decoupled the charged coulomb interactions by making the charges 0 in the [atoms] directive for the organic molecule. Then I ran energy minimization, equilibration (200 ps) and production run for 1 ns. I have varied the Lambda-vdw values :

free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
couple-moltype = Merocyanine ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
; 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
vdw_lambdas = 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.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
sc-sigma = 0.3
nstdhdl = 10

After choosing 21 Lambda values and running minimization, equilibration and production MD for all these Lambda values, I used gmx bar to calculate the solvation free energy. The delta G value for decoupling the organic molecule obtained from the above simulations is: -0.10 +/- 0.68 kJ/mol.

Is this value acceptable or it is erroreneous as the error bar is higher than the value itself. Also, how do I check if the free energy calculation is converged or do I need more Lambda values or do I need to run the simulations for longer time?
Any help would be appreciated, thank you

2

It looks like you need to read up on basic physics. You need to decide what accuracy you want. I would consider 0.7 kJ/mol a high accuracy in this context. You don’t write how long you simulated, but an error of 0.7 kJ/mol seems very small for such a large molecule.

Why did you manual zero the charges in the topology file? The couple option already zeros those for you.

3

Thank you for your reply. I was wondering that the free energy value is 0.1 kJ/mol whereas the error is six times larger than the free energy value and I was hoping if that is fine. Regarding the coulomb interactions, I was trying to follow a tutorial where the charges for the molecule is kept zero in the topology file and the mdp settings were provided accordingly. I might have to read more on the setting the mdp options.

4

I have run the equilibrium simulations for 200 ps and production MD for 1 ns for a little over 3000 atoms in the system for each lambda values.