Solvation Free Energy of NaCl

abhishekchemistry · January 20, 2023, 9:08am

GROMACS version: 2019.6
GROMACS modification: Yes/No
Here post your question

Hello, GROMACS users!

I have three questions regarding the solvation free energy calculation of NaCl using thermodynamic integration method. I have two systems: 0.01 M NaCl solution (1 NaCl in 5555 water) and 2.0 M NaCl (20 NaCl in 555 water). I am trying to calculate the solvation free energy of NaCl using thermodynamic integration and gmx bar method. I want to check the concentration dependance of solvation free energy. Here are my questions:

For both cases gmx bar is showing more or less same value for solvation free energy which is completely wrong according to the published data. I have no clue why this is happening.
Is it possible to decouple the interaction energy of one NaCl by just defining two dummy molecules (say DUMMY1 and DUMMY2 and Na belongs to DUMMY1 and Cl belongs to DUMMY2) and add them to “couple-moltype= DUMMY1, DUMMY2” to perform thermodynamic intregation in GROMACS?
I am a bit confused about “couple-intramol=yes/no” option here. Which one should I use?

I have pasted all the mdp, gro and top file here.

**## 0.01M NaCl files (glimpses): **

grompp.mdp

; Free energy control stuff
free_energy = yes
init_lambda_state = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
couple-moltype = DUMMY ; name of moleculetype to decouple
couple-lambda0 = vdw-q ; van der Waals and coulombic interactions
couple-lambda1 = none ; turn off everything
couple-intramol = yes

; 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 31 32 33 34 35 36 37 38 39 40
coul_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 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
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.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

; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; no: linear interpolation of Coulomb, yes: softcore Coulomb
sc-power = 1
sc-sigma = 0.3
nstdhdl = 1000

#conf.gro

5555SOL HW116664 0.402 1.434 5.069 -1.5551 -0.4714 -0.3539
5555SOL HW216665 0.511 1.349 4.982 0.4647 1.1323 0.5790
5556DUMMY NA_A16666 4.718 4.737 3.777 -0.4232 0.0243 -0.2703
5556DUMMY CL_A16667 2.193 2.749 0.481 -0.2597 0.0810 0.0706
5.50705 5.50705 5.50705

# topol.top:

[ defaults ]
;nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 no 0.0 0.0

;Include forcefield parameters
#include “./jc_oplsaa.ff/forcefield.itp”

[ atomtypes ]
;name mass charge ptype sigma epsilon
NA_A 22.98977 1.00000 A 0.21600000 1.47500000 ;dummy Na
CL_A 35.453 -1.00000 A 0.48300000 0.05349200 ;dummy Cl

[ nonbond_params ]
;i j func sigma epsilon

[ moleculetype ]
;molname nrexcl
DUMMY 2

[ atoms ]
;nr type resnr residue atom cgnr
1 NA_A 1 DUM NA_A 1
2 CL_A 1 DUM CL_A 1

; Include topology for ions
#include “./jc_oplsaa.ff/ions.itp”

; Include water topology
#include “./jc_oplsaa.ff/spce.itp”

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif

[ system ]
Water555-NaCl20

[ molecules ]
SOL 5555
DUMMY 1

**## 2.0M NaCl files (glimpses): **

# grompp.mdp:
same for 0.01M NaCl

#conf.gro
…
555SOL HW1 1664 1.474 0.445 1.930 0.7738 -0.3329 -0.5133
555SOL HW2 1665 1.312 0.423 1.925 0.6031 0.8510 -0.3764
556NA NA 1666 1.413 1.231 2.310 -0.0947 -0.4056 -0.2063
557NA NA 1667 0.200 1.276 1.347 0.0165 -0.0163 0.0772
558NA NA 1668 1.445 0.263 2.172 0.5565 -0.2146 0.3465
559NA NA 1669 0.975 0.125 0.624 -0.4964 -0.2147 -0.0458
560NA NA 1670 1.143 1.788 0.962 -0.4533 -0.1442 -0.0811
561NA NA 1671 1.260 1.900 1.236 -0.0278 -0.1462 0.3696
562NA NA 1672 1.901 0.690 1.547 0.6126 -0.1148 0.0389
563NA NA 1673 1.845 2.161 2.040 -0.2882 0.2941 -0.0610
564NA NA 1674 1.057 0.881 0.937 0.0523 -0.2782 -0.1406
565NA NA 1675 2.374 2.035 2.435 0.0308 -0.0964 -0.1752
566NA NA 1676 1.688 0.017 2.248 -0.1554 -0.1625 -0.1218
567NA NA 1677 1.756 0.435 0.434 0.3080 0.2392 0.1733
568NA NA 1678 1.323 2.287 2.387 -0.5513 0.4337 -0.2610
569NA NA 1679 2.068 0.245 1.962 -0.0904 0.1738 -0.0407
570NA NA 1680 1.962 1.183 0.611 -0.0532 -0.3497 -0.3122
571NA NA 1681 2.571 1.843 0.406 0.0291 -0.0419 0.5745
572NA NA 1682 0.475 0.112 1.676 0.0963 0.3003 -0.4348
573NA NA 1683 2.296 1.807 1.176 0.4238 0.0810 -0.2627
574NA NA 1684 0.266 1.834 2.106 0.1954 0.1261 0.0140
575NA NA 1685 0.968 0.990 2.141 0.1855 -0.4380 -0.0803
576CL CL 1687 1.773 1.087 0.486 -0.3116 0.4077 0.1034
577CL CL 1688 1.090 0.747 0.725 -0.2556 -0.0858 -0.2714
578CL CL 1689 1.611 1.500 2.033 -0.0733 0.1332 0.2154
579CL CL 1690 1.932 1.395 0.126 -0.0659 -0.5642 0.2723
580CL CL 1691 0.682 1.528 0.576 0.1249 -0.0542 0.0680
581CL CL 1692 0.766 0.544 1.998 -0.0538 -0.0371 0.2025
582CL CL 1693 0.185 1.472 1.538 -0.2161 -0.3380 -0.7710
583CL CL 1694 2.417 1.116 1.135 0.1265 0.3586 -0.2300
584CL CL 1695 0.310 1.638 1.958 -0.1292 0.3156 -0.1641
585CL CL 1696 1.861 0.385 1.896 0.2516 -0.0816 0.0174
586CL CL 1697 0.005 2.583 1.769 0.2247 0.3245 0.1556
587CL CL 1698 1.214 1.308 1.840 0.0745 0.0944 -0.0963
588CL CL 1699 1.125 1.644 0.444 -0.1223 0.0955 0.0590
589CL CL 1700 0.837 2.507 0.568 0.1852 0.4097 0.0244
590CL CL 1701 1.079 1.707 1.214 0.0614 -0.2608 0.2629
591CL CL 1702 1.511 0.116 2.403 -0.4012 -0.6571 -0.0118
592CL CL 1703 0.326 2.075 2.128 0.1719 -0.0552 0.2096
593CL CL 1704 1.756 2.355 0.080 0.0855 -0.1551 0.0488
594CL CL 1705 1.631 2.414 1.671 0.2301 0.0258 -0.1349
595CL CL 1706 2.146 1.445 0.933 -0.4356 -0.1190 -0.1822
596DUM NA_A 1686 1.898 0.862 2.114 0.1962 0.0222 0.3801
596DUM CL_A 1707 0.047 0.892 2.113 0.0672 0.0278 -0.0299
2.58723 2.58723 2.58723