Variance in chemical potential computation

dpo854 commented 3 years ago

Dear All, I noticed a peculiar behavior during Widom insertion for excess chemical potential (mu excess) computation of a salt pair. This peculiarity is in the variance in mu excess. Ideally as one adds more insertion cycles to the program (a cycle consists of a few reaction moves followed by a few MD runs. The latter is to rule out a correlation between the base configurations of successive insertions), the variance should either go down or stabilize to a small value. Here, I am seeing an abrupt increase in it after a few hundred cycles. The original variance was around 1.3 percent while the final is around 3 percent. The behavior is illustrated below in the output of the program. The system has two pairs: the first, cation=+3 and anion=-1, and the second, H+ and OH-. The jump is highlighted in bold.

Thanks, Debadutta Prusty

output (click to unroll)

mu_ex_pair 1 (-5.7285, +/- 0.8462) mu_ex_pair 2 (-1.3334, +/- 0.3292) mu_ex_pair 1 (-8.0667, +/- 0.9013) mu_ex_pair 2 (-1.2788, +/- 0.1792) mu_ex_pair 1 (-7.4791, +/- 0.6436) mu_ex_pair 2 (-1.3537, +/- 0.1139) mu_ex_pair 1 (-7.3687, +/- 0.4984) mu_ex_pair 2 (-1.3125, +/- 0.0838) mu_ex_pair 1 (-7.4090, +/- 0.4011) mu_ex_pair 2 (-1.4156, +/- 0.0966) mu_ex_pair 1 (-7.2962, +/- 0.3452) mu_ex_pair 2 (-1.3608, +/- 0.0804) mu_ex_pair 1 (-7.4538, +/- 0.3172) mu_ex_pair 2 (-1.4192, +/- 0.0782) mu_ex_pair 1 (-7.2861, +/- 0.2922) mu_ex_pair 2 (-1.4086, +/- 0.0691) mu_ex_pair 1 (-7.3382, +/- 0.2628) mu_ex_pair 2 (-1.3987, +/- 0.0627) mu_ex_pair 1 (-7.2547, +/- 0.2435) mu_ex_pair 2 (-1.3859, +/- 0.0579) mu_ex_pair 1 (-7.2012, +/- 0.2259) mu_ex_pair 2 (-1.3676, +/- 0.0533) mu_ex_pair 1 (-7.1692, +/- 0.2123) mu_ex_pair 2 (-1.3966, +/- 0.0577) mu_ex_pair 1 (-7.1177, +/- 0.2000) mu_ex_pair 2 (-1.3731, +/- 0.0539) mu_ex_pair 1 (-7.0849, +/- 0.1889) mu_ex_pair 2 (-1.3720, +/- 0.0530) mu_ex_pair 1 (-7.1082, +/- 0.1860) mu_ex_pair 2 (-1.3758, +/- 0.0512) mu_ex_pair 1 (-7.0943, +/- 0.1772) mu_ex_pair 2 (-1.3660, +/- 0.0487) mu_ex_pair 1 (-7.1035, +/- 0.1693) mu_ex_pair 2 (-1.3479, +/- 0.0462) mu_ex_pair 1 (-7.0838, +/- 0.1621) mu_ex_pair 2 (-1.3504, +/- 0.0440) mu_ex_pair 1 (-7.0401, +/- 0.1565) mu_ex_pair 2 (-1.3589, +/- 0.0429) mu_ex_pair 1 (-6.9861, +/- 0.1520) mu_ex_pair 2 (-1.3664, +/- 0.0421) mu_ex_pair 1 (-7.0123, +/- 0.1476) mu_ex_pair 2 (-1.3628, +/- 0.0408) mu_ex_pair 1 (-7.0181, +/- 0.1429) mu_ex_pair 2 (-1.3878, +/- 0.0416) mu_ex_pair 1 (-7.0443, +/- 0.1428) mu_ex_pair 2 (-1.3965, +/- 0.0435) mu_ex_pair 1 (-7.0328, +/- 0.1383) mu_ex_pair 2 (-1.3922, +/- 0.0421) mu_ex_pair 1 (-7.0136, +/- 0.1354) mu_ex_pair 2 (-1.3884, +/- 0.0408) mu_ex_pair 1 (-7.0296, +/- 0.1389) mu_ex_pair 2 (-1.3769, +/- 0.0395) mu_ex_pair 1 (-7.0533, +/- 0.1347) mu_ex_pair 2 (-1.3685, +/- 0.0384) mu_ex_pair 1 (-7.0667, +/- 0.1333) mu_ex_pair 2 (-1.3697, +/- 0.0374) mu_ex_pair 1 (-7.0997, +/- 0.1320) mu_ex_pair 2 (-1.3871, +/- 0.0379) mu_ex_pair 1 (-7.1371, +/- 0.1316) mu_ex_pair 2 (-1.3780, +/- 0.0368) mu_ex_pair 1 (-7.1672, +/- 0.1302) mu_ex_pair 2 (-1.3794, +/- 0.0361) mu_ex_pair 1 (-7.1885, +/- 0.1270) mu_ex_pair 2 (-1.3783, +/- 0.0352) mu_ex_pair 1 (-7.1814, +/- 0.1248) mu_ex_pair 2 (-1.3750, +/- 0.0344) mu_ex_pair 1 (-7.1745, +/- 0.1221) mu_ex_pair 2 (-1.3886, +/- 0.0385) mu_ex_pair 1 (-7.1500, +/- 0.1198) mu_ex_pair 2 (-1.3976, +/- 0.0378) mu_ex_pair 1 (-7.1284, +/- 0.1176) mu_ex_pair 2 (-1.3919, +/- 0.0369) mu_ex_pair 1 (-7.1371, +/- 0.1166) mu_ex_pair 2 (-1.4190, +/- 0.0385) mu_ex_pair 1 (-7.1763, +/- 0.1221) mu_ex_pair 2 (-1.4180, +/- 0.0385) mu_ex_pair 1 (-7.1645, +/- 0.1198) mu_ex_pair 2 (-1.4149, +/- 0.0377) mu_ex_pair 1 (-7.1456, +/- 0.1178) mu_ex_pair 2 (-1.4075, +/- 0.0368) mu_ex_pair 1 (-7.1197, +/- 0.1162) mu_ex_pair 2 (-1.4064, +/- 0.0361) mu_ex_pair 1 (-7.1162, +/- 0.1166) mu_ex_pair 2 (-1.4079, +/- 0.0354) mu_ex_pair 1 (-7.1085, +/- 0.1148) mu_ex_pair 2 (-1.4057, +/- 0.0351) mu_ex_pair 1 (-7.2933, +/- 0.2101) mu_ex_pair 2 (-1.4018, +/- 0.0344) mu_ex_pair 1 (-7.3202, +/- 0.2051) mu_ex_pair 2 (-1.4078, +/- 0.0341) mu_ex_pair 1 (-7.3022, +/- 0.2020) mu_ex_pair 2 (-1.4090, +/- 0.0337) mu_ex_pair 1 (-7.2906, +/- 0.1988) mu_ex_pair 2 (-1.4120, +/- 0.0332) mu_ex_pair 1 (-7.2936, +/- 0.1954) mu_ex_pair 2 (-1.4131, +/- 0.0326) mu_ex_pair 1 (-7.2970, +/- 0.1913) mu_ex_pair 2 (-1.4081, +/- 0.0321) mu_ex_pair 1 (-7.3089, +/- 0.1878) mu_ex_pair 2 (-1.4066, +/- 0.0315) mu_ex_pair 1 (-7.2956, +/- 0.1853) mu_ex_pair 2 (-1.4046, +/- 0.0310) mu_ex_pair 1 (-7.2867, +/- 0.1825) mu_ex_pair 2 (-1.4101, +/- 0.0309) mu_ex_pair 1 (-7.2739, +/- 0.1801) mu_ex_pair 2 (-1.4141, +/- 0.0307) mu_ex_pair 1 (-7.2802, +/- 0.1775) mu_ex_pair 2 (-1.4156, +/- 0.0304) mu_ex_pair 1 (-7.2740, +/- 0.1749) mu_ex_pair 2 (-1.4199, +/- 0.0303) mu_ex_pair 1 (-7.2576, +/- 0.1729) mu_ex_pair 2 (-1.4175, +/- 0.0299) mu_ex_pair 1 (-7.2590, +/- 0.1700) mu_ex_pair 2 (-1.4131, +/- 0.0295) mu_ex_pair 1 (-7.2537, +/- 0.1677) mu_ex_pair 2 (-1.4084, +/- 0.0290) mu_ex_pair 1 (-7.3195, +/- 0.1768) mu_ex_pair 2 (-1.4035, +/- 0.0287) mu_ex_pair 1 (-7.3187, +/- 0.1746) mu_ex_pair 2 (-1.4005, +/- 0.0283) mu_ex_pair 1 (-7.3334, +/- 0.1727) mu_ex_pair 2 (-1.3966, +/- 0.0279) mu_ex_pair 1 (-7.3199, +/- 0.1709) mu_ex_pair 2 (-1.3919, +/- 0.0275) mu_ex_pair 1 (-7.3478, +/- 0.1695) mu_ex_pair 2 (-1.3922, +/- 0.0273) mu_ex_pair 1 (-7.3345, +/- 0.1677) mu_ex_pair 2 (-1.3961, +/- 0.0275) mu_ex_pair 1 (-7.3375, +/- 0.1653) mu_ex_pair 2 (-1.3989, +/- 0.0272) mu_ex_pair 1 (-7.3310, +/- 0.1634) mu_ex_pair 2 (-1.3969, +/- 0.0269) mu_ex_pair 1 (-7.3600, +/- 0.1638) mu_ex_pair 2 (-1.4055, +/- 0.0272) mu_ex_pair 1 (-7.3551, +/- 0.1619) mu_ex_pair 2 (-1.4084, +/- 0.0275) mu_ex_pair 1 (-7.3383, +/- 0.1606) mu_ex_pair 2 (-1.4068, +/- 0.0272) mu_ex_pair 1 (-7.3197, +/- 0.1595) mu_ex_pair 2 (-1.4050, +/- 0.0269) mu_ex_pair 1 (-7.3092, +/- 0.1579) mu_ex_pair 2 (-1.4014, +/- 0.0266) mu_ex_pair 1 (-7.3064, +/- 0.1560) mu_ex_pair 2 (-1.4014, +/- 0.0266) mu_ex_pair 1 (-7.3130, +/- 0.1536) mu_ex_pair 2 (-1.4031, +/- 0.0265) mu_ex_pair 1 (-7.2999, +/- 0.1524) mu_ex_pair 2 (-1.4001, +/- 0.0262) mu_ex_pair 1 (-7.3211, +/- 0.1518) mu_ex_pair 2 (-1.3959, +/- 0.0259) mu_ex_pair 1 (-7.3226, +/- 0.1500) mu_ex_pair 2 (-1.3958, +/- 0.0256) mu_ex_pair 1 (-7.3136, +/- 0.1486) mu_ex_pair 2 (-1.3923, +/- 0.0254) mu_ex_pair 1 (-7.3179, +/- 0.1466) mu_ex_pair 2 (-1.3904, +/- 0.0252) mu_ex_pair 1 (-7.3082, +/- 0.1453) mu_ex_pair 2 (-1.3887, +/- 0.0249) mu_ex_pair 1 (-7.3021, +/- 0.1439) mu_ex_pair 2 (-1.3871, +/- 0.0247) mu_ex_pair 1 (-7.3072, +/- 0.1424) mu_ex_pair 2 (-1.3870, +/- 0.0244) mu_ex_pair 1 (-7.3024, +/- 0.1410) mu_ex_pair 2 (-1.3882, +/- 0.0243) mu_ex_pair 1 (-7.3110, +/- 0.1391) mu_ex_pair 2 (-1.3847, +/- 0.0240) mu_ex_pair 1 (-7.3034, +/- 0.1379) mu_ex_pair 2 (-1.3838, +/- 0.0238) mu_ex_pair 1 (-7.2928, +/- 0.1369) mu_ex_pair 2 (-1.3868, +/- 0.0239) mu_ex_pair 1 (-7.2982, +/- 0.1354) mu_ex_pair 2 (-1.3859, +/- 0.0237) mu_ex_pair 1 (-7.2976, +/- 0.1340) mu_ex_pair 2 (-1.3824, +/- 0.0235) mu_ex_pair 1 (-7.2914, +/- 0.1329) mu_ex_pair 2 (-1.3811, +/- 0.0233) mu_ex_pair 1 (-7.3377, +/- 0.1375) mu_ex_pair 2 (-1.3849, +/- 0.0231) mu_ex_pair 1 (-7.3320, +/- 0.1363) mu_ex_pair 2 (-1.3847, +/- 0.0229) mu_ex_pair 1 (-7.3493, +/- 0.1358) mu_ex_pair 2 (-1.3882, +/- 0.0229) mu_ex_pair 1 (-7.3573, +/- 0.1344) mu_ex_pair 2 (-1.3890, +/- 0.0227) mu_ex_pair 1 (-7.3682, +/- 0.1331) mu_ex_pair 2 (-1.3865, +/- 0.0226) mu_ex_pair 1 (-7.3594, +/- 0.1322) mu_ex_pair 2 (-1.3841, +/- 0.0224) mu_ex_pair 1 (-7.3539, +/- 0.1311) mu_ex_pair 2 (-1.3834, +/- 0.0222) mu_ex_pair 1 (-7.3500, +/- 0.1300) mu_ex_pair 2 (-1.3883, +/- 0.0226) mu_ex_pair 1 (-7.3457, +/- 0.1290) mu_ex_pair 2 (-1.3883, +/- 0.0224) mu_ex_pair 1 (-7.3681, +/- 0.1282) mu_ex_pair 2 (-1.3892, +/- 0.0222) mu_ex_pair 1 (-7.3629, +/- 0.1273) mu_ex_pair 2 (-1.3905, +/- 0.0222) mu_ex_pair 1 (-7.3530, +/- 0.1265) mu_ex_pair 2 (-1.3909, +/- 0.0220) mu_ex_pair 1 (-7.3502, +/- 0.1255) mu_ex_pair 2 (-1.3917, +/- 0.0219) mu_ex_pair 1 (-7.3489, +/- 0.1244) mu_ex_pair 2 (-1.3930, +/- 0.0217) mu_ex_pair 1 (-7.3502, +/- 0.1236) mu_ex_pair 2 (-1.3940, +/- 0.0218) mu_ex_pair 1 (-7.3417, +/- 0.1228) mu_ex_pair 2 (-1.3925, +/- 0.0216) mu_ex_pair 1 (-7.3334, +/- 0.1221) mu_ex_pair 2 (-1.3905, +/- 0.0215) mu_ex_pair 1 (-7.3340, +/- 0.1210) mu_ex_pair 2 (-1.3883, +/- 0.0213) mu_ex_pair 1 (-7.3287, +/- 0.1202) mu_ex_pair 2 (-1.3897, +/- 0.0212) mu_ex_pair 1 (-7.3207, +/- 0.1195) mu_ex_pair 2 (-1.3881, +/- 0.0211) mu_ex_pair 1 (-7.3116, +/- 0.1189) mu_ex_pair 2 (-1.3907, +/- 0.0210) mu_ex_pair 1 (-7.3054, +/- 0.1181) mu_ex_pair 2 (-1.3905, +/- 0.0210) mu_ex_pair 1 (-7.3129, +/- 0.1174) mu_ex_pair 2 (-1.3875, +/- 0.0208) mu_ex_pair 1 (-7.3102, +/- 0.1165) mu_ex_pair 2 (-1.3852, +/- 0.0207) mu_ex_pair 1 (-7.3025, +/- 0.1159) mu_ex_pair 2 (-1.3855, +/- 0.0205) mu_ex_pair 1 (-7.3020, +/- 0.1150) mu_ex_pair 2 (-1.3884, +/- 0.0205) mu_ex_pair 1 (-7.3144, +/- 0.1142) mu_ex_pair 2 (-1.3853, +/- 0.0204) mu_ex_pair 1 (-7.3082, +/- 0.1135) mu_ex_pair 2 (-1.3870, +/- 0.0203) mu_ex_pair 1 (-7.2994, +/- 0.1129) mu_ex_pair 2 (-1.3848, +/- 0.0201) mu_ex_pair 1 (-7.3002, +/- 0.1121) mu_ex_pair 2 (-1.3859, +/- 0.0201) mu_ex_pair 1 (-7.2985, +/- 0.1114) mu_ex_pair 2 (-1.3851, +/- 0.0199) mu_ex_pair 1 (-7.3251, +/- 0.1140) mu_ex_pair 2 (-1.3832, +/- 0.0198) mu_ex_pair 1 (-7.3223, +/- 0.1133) mu_ex_pair 2 (-1.3862, +/- 0.0198) mu_ex_pair 1 (-7.3178, +/- 0.1126) mu_ex_pair 2 (-1.3904, +/- 0.0198) mu_ex_pair 1 (-7.3118, +/- 0.1120) mu_ex_pair 2 (-1.3895, +/- 0.0197) mu_ex_pair 1 (-7.3076, +/- 0.1113) mu_ex_pair 2 (-1.3915, +/- 0.0197) mu_ex_pair 1 (-7.3021, +/- 0.1107) mu_ex_pair 2 (-1.3923, +/- 0.0196) mu_ex_pair 1 (-7.3196, +/- 0.1116) mu_ex_pair 2 (-1.3920, +/- 0.0194) mu_ex_pair 1 (-7.3449, +/- 0.1131) mu_ex_pair 2 (-1.3922, +/- 0.0194) mu_ex_pair 1 (-7.3399, +/- 0.1125) mu_ex_pair 2 (-1.3912, +/- 0.0192) mu_ex_pair 1 (-7.3297, +/- 0.1121) mu_ex_pair 2 (-1.3944, +/- 0.0193) mu_ex_pair 1 (-7.3247, +/- 0.1115) mu_ex_pair 2 (-1.3936, +/- 0.0192) mu_ex_pair 1 (-7.3265, +/- 0.1107) mu_ex_pair 2 (-1.3922, +/- 0.0191) mu_ex_pair 1 (-7.3271, +/- 0.1100) mu_ex_pair 2 (-1.3909, +/- 0.0190) mu_ex_pair 1 (-7.3338, +/- 0.1096) mu_ex_pair 2 (-1.3888, +/- 0.0188) mu_ex_pair 1 (-7.3685, +/- 0.1115) mu_ex_pair 2 (-1.3908, +/- 0.0188) mu_ex_pair 1 (-7.3608, +/- 0.1110) mu_ex_pair 2 (-1.3906, +/- 0.0187) mu_ex_pair 1 (-7.3570, +/- 0.1104) mu_ex_pair 2 (-1.3890, +/- 0.0186) mu_ex_pair 1 (-7.3530, +/- 0.1098) mu_ex_pair 2 (-1.3867, +/- 0.0185) mu_ex_pair 1 (-7.3646, +/- 0.1099) mu_ex_pair 2 (-1.3900, +/- 0.0185) mu_ex_pair 1 (-7.3609, +/- 0.1093) mu_ex_pair 2 (-1.3900, +/- 0.0184) mu_ex_pair 1 (-7.3620, +/- 0.1087) mu_ex_pair 2 (-1.3913, +/- 0.0185) mu_ex_pair 1 (-7.3557, +/- 0.1082) mu_ex_pair 2 (-1.3898, +/- 0.0184) mu_ex_pair 1 (-7.3493, +/- 0.1077) mu_ex_pair 2 (-1.3890, +/- 0.0183) mu_ex_pair 1 (-7.3483, +/- 0.1070) mu_ex_pair 2 (-1.3868, +/- 0.0181) mu_ex_pair 1 (-7.3427, +/- 0.1066) mu_ex_pair 2 (-1.3883, +/- 0.0183) mu_ex_pair 1 (-7.3404, +/- 0.1060) mu_ex_pair 2 (-1.3864, +/- 0.0181) mu_ex_pair 1 (-7.3381, +/- 0.1054) mu_ex_pair 2 (-1.3851, +/- 0.0181) mu_ex_pair 1 (-7.3397, +/- 0.1048) mu_ex_pair 2 (-1.3857, +/- 0.0180) mu_ex_pair 1 (-7.3323, +/- 0.1044) mu_ex_pair 2 (-1.3868, +/- 0.0180) mu_ex_pair 1 (-7.3293, +/- 0.1039) mu_ex_pair 2 (-1.3854, +/- 0.0179) mu_ex_pair 1 (-7.3241, +/- 0.1034) mu_ex_pair 2 (-1.3890, +/- 0.0180) mu_ex_pair 1 (-7.3231, +/- 0.1028) mu_ex_pair 2 (-1.3879, +/- 0.0179) mu_ex_pair 1 (-7.3280, +/- 0.1021) mu_ex_pair 2 (-1.3859, +/- 0.0178) mu_ex_pair 1 (-7.3302, +/- 0.1017) mu_ex_pair 2 (-1.3846, +/- 0.0177) mu_ex_pair 1 (-7.3290, +/- 0.1011) mu_ex_pair 2 (-1.3878, +/- 0.0178) mu_ex_pair 1 (-7.3251, +/- 0.1007) mu_ex_pair 2 (-1.3871, +/- 0.0177) mu_ex_pair 1 (-7.3186, +/- 0.1003) mu_ex_pair 2 (-1.3873, +/- 0.0176) mu_ex_pair 1 (-7.3097, +/- 0.1000) mu_ex_pair 2 (-1.3857, +/- 0.0175) mu_ex_pair 1 (-7.3046, +/- 0.0996) mu_ex_pair 2 (-1.3858, +/- 0.0174) mu_ex_pair 1 (-7.3101, +/- 0.0990) mu_ex_pair 2 (-1.3848, +/- 0.0173) mu_ex_pair 1 (-7.3122, +/- 0.0984) mu_ex_pair 2 (-1.3847, +/- 0.0173) mu_ex_pair 1 (-7.3102, +/- 0.0979) mu_ex_pair 2 (-1.3838, +/- 0.0172) mu_ex_pair 1 (-7.3069, +/- 0.0975) mu_ex_pair 2 (-1.3827, +/- 0.0171) mu_ex_pair 1 (-7.3054, +/- 0.0970) mu_ex_pair 2 (-1.3816, +/- 0.0170) mu_ex_pair 1 (-7.3049, +/- 0.0965) mu_ex_pair 2 (-1.3824, +/- 0.0170) mu_ex_pair 1 (-7.3008, +/- 0.0961) mu_ex_pair 2 (-1.3804, +/- 0.0169) mu_ex_pair 1 (-7.3010, +/- 0.0956) mu_ex_pair 2 (-1.3803, +/- 0.0169) mu_ex_pair 1 (-7.3046, +/- 0.0950) mu_ex_pair 2 (-1.3794, +/- 0.0168) mu_ex_pair 1 (-7.2963, +/- 0.0947) mu_ex_pair 2 (-1.3791, +/- 0.0167) mu_ex_pair 1 (-7.2911, +/- 0.0944) mu_ex_pair 2 (-1.3773, +/- 0.0166) mu_ex_pair 1 (-7.2918, +/- 0.0938) mu_ex_pair 2 (-1.3767, +/- 0.0165) mu_ex_pair 1 (-7.2837, +/- 0.0936) mu_ex_pair 2 (-1.3780, +/- 0.0167) mu_ex_pair 1 (-7.2832, +/- 0.0931) mu_ex_pair 2 (-1.3781, +/- 0.0166) mu_ex_pair 1 (-7.2992, +/- 0.0939) mu_ex_pair 2 (-1.3800, +/- 0.0166) mu_ex_pair 1 (-7.3013, +/- 0.0934) mu_ex_pair 2 (-1.3813, +/- 0.0165) mu_ex_pair 1 (-7.3014, +/- 0.0929) mu_ex_pair 2 (-1.3803, +/- 0.0164) mu_ex_pair 1 (-7.2979, +/- 0.0926) mu_ex_pair 2 (-1.3779, +/- 0.0164) mu_ex_pair 1 (-7.3109, +/- 0.0927) mu_ex_pair 2 (-1.3769, +/- 0.0163) mu_ex_pair 1 (-7.3074, +/- 0.0923) mu_ex_pair 2 (-1.3762, +/- 0.0162) mu_ex_pair 1 (-7.3024, +/- 0.0920) mu_ex_pair 2 (-1.3751, +/- 0.0162) **mu_ex_pair 1 (-7.3001, +/- 0.0916) mu_ex_pair 2 (-1.3755, +/- 0.0161) mu_ex_pair 1 (-7.7106, +/- 0.3908) mu_ex_pair 2 (-1.3750, +/- 0.0160)** mu_ex_pair 1 (-7.7121, +/- 0.3885) mu_ex_pair 2 (-1.3744, +/- 0.0160) mu_ex_pair 1 (-7.7080, +/- 0.3870) mu_ex_pair 2 (-1.3772, +/- 0.0163) mu_ex_pair 1 (-7.7024, +/- 0.3858) mu_ex_pair 2 (-1.3760, +/- 0.0162) mu_ex_pair 1 (-7.6972, +/- 0.3845) mu_ex_pair 2 (-1.3757, +/- 0.0162) mu_ex_pair 1 (-7.6935, +/- 0.3830) mu_ex_pair 2 (-1.3751, +/- 0.0161) mu_ex_pair 1 (-7.6918, +/- 0.3812) mu_ex_pair 2 (-1.3763, +/- 0.0160) mu_ex_pair 1 (-7.6949, +/- 0.3788) mu_ex_pair 2 (-1.3750, +/- 0.0160) mu_ex_pair 1 (-7.6895, +/- 0.3776) mu_ex_pair 2 (-1.3768, +/- 0.0160) mu_ex_pair 1 (-7.7116, +/- 0.3732) mu_ex_pair 2 (-1.3773, +/- 0.0159) mu_ex_pair 1 (-7.7050, +/- 0.3723) mu_ex_pair 2 (-1.3786, +/- 0.0159) mu_ex_pair 1 (-7.7023, +/- 0.3707) mu_ex_pair 2 (-1.3788, +/- 0.0158) mu_ex_pair 1 (-7.6961, +/- 0.3697) mu_ex_pair 2 (-1.3786, +/- 0.0158) mu_ex_pair 1 (-7.7216, +/- 0.3650) mu_ex_pair 2 (-1.3787, +/- 0.0157) mu_ex_pair 1 (-7.7170, +/- 0.3638) mu_ex_pair 2 (-1.3788, +/- 0.0157) mu_ex_pair 1 (-7.7167, +/- 0.3620) mu_ex_pair 2 (-1.3785, +/- 0.0156) mu_ex_pair 1 (-7.7122, +/- 0.3608) mu_ex_pair 2 (-1.3787, +/- 0.0156) mu_ex_pair 1 (-7.7114, +/- 0.3592) mu_ex_pair 2 (-1.3778, +/- 0.0155) mu_ex_pair 1 (-7.7227, +/- 0.3560) mu_ex_pair 2 (-1.3822, +/- 0.0156) mu_ex_pair 1 (-7.7605, +/- 0.3510) mu_ex_pair 2 (-1.3828, +/- 0.0155) mu_ex_pair 1 (-7.7560, +/- 0.3498) mu_ex_pair 2 (-1.3834, +/- 0.0155) mu_ex_pair 1 (-7.7505, +/- 0.3489) mu_ex_pair 2 (-1.3832, +/- 0.0154) mu_ex_pair 1 (-7.7466, +/- 0.3477) mu_ex_pair 2 (-1.3825, +/- 0.0154) mu_ex_pair 1 (-7.7400, +/- 0.3469) mu_ex_pair 2 (-1.3811, +/- 0.0153) mu_ex_pair 1 (-7.7382, +/- 0.3455) mu_ex_pair 2 (-1.3799, +/- 0.0153) mu_ex_pair 1 (-7.7373, +/- 0.3439) mu_ex_pair 2 (-1.3788, +/- 0.0152) mu_ex_pair 1 (-7.7381, +/- 0.3422) mu_ex_pair 2 (-1.3781, +/- 0.0151) mu_ex_pair 1 (-7.7357, +/- 0.3409) mu_ex_pair 2 (-1.3776, +/- 0.0151) mu_ex_pair 1 (-7.7333, +/- 0.3396) mu_ex_pair 2 (-1.3799, +/- 0.0151) mu_ex_pair 1 (-7.7272, +/- 0.3388) mu_ex_pair 2 (-1.3786, +/- 0.0151) mu_ex_pair 1 (-7.7236, +/- 0.3377) mu_ex_pair 2 (-1.3803, +/- 0.0152) mu_ex_pair 1 (-7.7202, +/- 0.3365) mu_ex_pair 2 (-1.3797, +/- 0.0152) mu_ex_pair 1 (-7.7185, +/- 0.3352) mu_ex_pair 2 (-1.3815, +/- 0.0153) mu_ex_pair 1 (-7.7142, +/- 0.3342) mu_ex_pair 2 (-1.3815, +/- 0.0152) mu_ex_pair 1 (-7.7198, +/- 0.3320) mu_ex_pair 2 (-1.3831, +/- 0.0152) mu_ex_pair 1 (-7.7172, +/- 0.3308) mu_ex_pair 2 (-1.3828, +/- 0.0152) mu_ex_pair 1 (-7.7137, +/- 0.3297) mu_ex_pair 2 (-1.3827, +/- 0.0151) mu_ex_pair 1 (-7.7097, +/- 0.3287) mu_ex_pair 2 (-1.3814, +/- 0.0151) mu_ex_pair 1 (-7.7065, +/- 0.3276) mu_ex_pair 2 (-1.3832, +/- 0.0153) mu_ex_pair 1 (-7.7129, +/- 0.3255) mu_ex_pair 2 (-1.3821, +/- 0.0152) mu_ex_pair 1 (-7.7080, +/- 0.3247) mu_ex_pair 2 (-1.3830, +/- 0.0152) mu_ex_pair 1 (-7.7077, +/- 0.3233) mu_ex_pair 2 (-1.3872, +/- 0.0157) mu_ex_pair 1 (-7.7051, +/- 0.3222) mu_ex_pair 2 (-1.3874, +/- 0.0156) mu_ex_pair 1 (-7.7003, +/- 0.3213) mu_ex_pair 2 (-1.3883, +/- 0.0156) mu_ex_pair 1 (-7.6975, +/- 0.3203) mu_ex_pair 2 (-1.3881, +/- 0.0155) mu_ex_pair 1 (-7.6930, +/- 0.3194) mu_ex_pair 2 (-1.3871, +/- 0.0155) mu_ex_pair 1 (-7.6910, +/- 0.3183) mu_ex_pair 2 (-1.3859, +/- 0.0154) mu_ex_pair 1 (-7.6862, +/- 0.3175) mu_ex_pair 2 (-1.3854, +/- 0.0154) mu_ex_pair 1 (-7.6884, +/- 0.3159) mu_ex_pair 2 (-1.3848, +/- 0.0153) mu_ex_pair 1 (-7.6846, +/- 0.3150) mu_ex_pair 2 (-1.3861, +/- 0.0154) mu_ex_pair 1 (-7.6804, +/- 0.3142) mu_ex_pair 2 (-1.3853, +/- 0.0153) mu_ex_pair 1 (-7.6784, +/- 0.3131) mu_ex_pair 2 (-1.3840, +/- 0.0153) mu_ex_pair 1 (-7.6765, +/- 0.3120) mu_ex_pair 2 (-1.3850, +/- 0.0153) mu_ex_pair 1 (-7.6711, +/- 0.3113) mu_ex_pair 2 (-1.3865, +/- 0.0152) mu_ex_pair 1 (-7.6713, +/- 0.3100) mu_ex_pair 2 (-1.3861, +/- 0.0152) mu_ex_pair 1 (-7.6672, +/- 0.3092) mu_ex_pair 2 (-1.3868, +/- 0.0152) mu_ex_pair 1 (-7.6622, +/- 0.3085) mu_ex_pair 2 (-1.3868, +/- 0.0151) mu_ex_pair 1 (-7.6579, +/- 0.3077) mu_ex_pair 2 (-1.3872, +/- 0.0151) mu_ex_pair 1 (-7.6542, +/- 0.3069) mu_ex_pair 2 (-1.3874, +/- 0.0151) mu_ex_pair 1 (-7.6506, +/- 0.3060) mu_ex_pair 2 (-1.3859, +/- 0.0150) mu_ex_pair 1 (-7.6478, +/- 0.3051) mu_ex_pair 2 (-1.3849, +/- 0.0150) mu_ex_pair 1 (-7.6450, +/- 0.3042) mu_ex_pair 2 (-1.3850, +/- 0.0149) mu_ex_pair 1 (-7.6416, +/- 0.3033) mu_ex_pair 2 (-1.3861, +/- 0.0149) mu_ex_pair 1 (-7.6386, +/- 0.3024) mu_ex_pair 2 (-1.3855, +/- 0.0149) mu_ex_pair 1 (-7.6347, +/- 0.3017) mu_ex_pair 2 (-1.3858, +/- 0.0148) mu_ex_pair 1 (-7.6340, +/- 0.3006) mu_ex_pair 2 (-1.3849, +/- 0.0148) mu_ex_pair 1 (-7.6294, +/- 0.2999) mu_ex_pair 2 (-1.3840, +/- 0.0147) mu_ex_pair 1 (-7.6359, +/- 0.2980) mu_ex_pair 2 (-1.3832, +/- 0.0147) mu_ex_pair 1 (-7.6388, +/- 0.2965) mu_ex_pair 2 (-1.3822, +/- 0.0146) mu_ex_pair 1 (-7.6339, +/- 0.2959) mu_ex_pair 2 (-1.3812, +/- 0.0146)

jonaslandsgesell commented 3 years ago

There are several thoughts I want to add:

1) Widom insertion how it is implemented, is only valid for the canonical ensemble (no particle fluctuations after the Widom insertion step is finished). Do not use the implemented Widom insertion method e.g. in the reaction ensemble/constant pH ensemble/ grand canonical ensemble (as it is also stated in the documentation) - one would need to use modifications like https://www.tandfonline.com/doi/abs/10.1080/00268978800100743

2) Widoms method performs temporary insertions at completely arbitrary places in the simulation box (brute force sampling). The jumps you see for the reported estimate of the chemical potential and its error are (to my experience) due to the fact that rarely an insertion of an ion happens very close to another ion. This configuration has a very different energy and therfore contributes a lot to the estimate of average change of energy due to an insertion of an ion (most other configurations result in a smaller energy change). Therefore, the variance of the potential energy change estimate jumps (as well as the chemical potential itself). This problem becomes more pronounced for multivalent ions which have even higher interaction energies.

3) Do only use the provided error of the chemical potential as a rough estimate of the error of the chemical potential estimate. It is computed using the standard error (https://en.m.wikipedia.org/wiki/Standard_error) of the recorded changes in potential energies (due to inserting an additional particle). Possible correlations are neglected (but should also not be too high in brute force sampling of configurations). Neglecting correlations can result in underestimating the real error, https://en.m.wikipedia.org/wiki/Markov_chain_central_limit_theorem

PS: the values which are returned by the Widom insertion method are

the excess chemical potentials of the inserted ion (pairs): mu_ex,I=-kT ln( <exp(-beta Delta E_pot)>). Note the average is computed inside the logarithm!
the error estimate of the excess chemical potential (computed via propagation of error using above formula and the standard error of the mean potential energy change ) -- this is not the variance of the chemical potential...

jonaslandsgesell commented 3 years ago

PPS: Sampling states where the inserted particle is by chance close to an already existing one is rare, for the basic problem see also: https://www.ling.upenn.edu/courses/cogs502/LNRE.html To overcome these problems you can maybe use/ develop an advanced free energy sampling method.

jonaslandsgesell commented 3 years ago

Actually you could maybe adapt the Wang-Landau reaction ensemble method (without educt types) to make grand-canocial moves and obtain excess free energy estimations for inserting 1,2,3,4,... Ion pairs (the number of inserted ion pairs M is the collective variable, as well as the potential energy as another collective variable) from which you can compute the excess chemical potential mu_excess = dF_excess/dN evaluated at (M=< N >,V,T) where you could compute < N>=c*V from some concentration c for which you are interested in. I think this should work - one could publish this idea I think. This could be an interesting research question

dpo854 commented 3 years ago

Unfortunately, the chemical potential computation is just a small part of the problem I am trying to address. It is used as an input to another grand canonical simulation. So, I am happy with the value the current Widom insertion produces if the error is not high without these rare events.

By the way, the simulation is indeed in the canonical ensemble with 400 particles.

jonaslandsgesell commented 3 years ago

For sure the classical widom insertion method is accurate enough for what you want to do. In order to be on the safe side, just run a grand canonical simulation with the excess chemical potential for a given conce traction you determined as input and some random initial number of particles in the box. Then check whether the concentrations you want to achieve are actually achieved in the simulation box after equilibration.

dpo854 commented 3 years ago

Thanks for the clarification. I have a question regarding the grand canonical simulation though. Is the chemical potential that we get from Widom's insertion functionality beta mu or just simply mu? This distinction is important for grand canonical simulations since the input to the reaction ensemble asks for the chemical potential normalized by the temperature. So, if the output from Widom is already normalized by beta and I further divide it by the temperature, the grand canonical simulation can be faulty. I ask that because the particle number is increasing in my grand canonical simulation indefinitely with time.

jonaslandsgesell commented 3 years ago

Hi Debadutta,

Is the chemical potential that we get from Widom's insertion functionality beta mu or just simply mu Widom's insertion method is a method to compute the excess part of the chemical potential and therefore has dimensions of energy.

You can see here: https://github.com/espressomd/espresso/blob/ab38289d308df96c74383d5212b0a4e51b392477/src/core/reaction_ensemble.cpp#L1614 that the value of the excess chemical potential is returned including the temperature factor (k_BT) as it is documented here: http://espressomd.org/html/doc/advanced_methods.html#widom-insertion-for-homogeneous-systems

This distinction is important for grand canonical simulations since the input to the reaction ensemble asks for the chemical potential normalized by the temperature

It is a bit more involved, I would say. The reaction ensemble can be "misused" to mimic a grand-canonical simulation. This is possible and totally fine because the reaction ensemble is derived from the grand canonical ensemble: https://www.tandfonline.com/doi/abs/10.1080/08927020801986564

The input to the reaction ensemble is:

particle types of the to be created/deleted particles
their stochiometry
the volume V
temperature (kT)
as well as the "reaction constant" gamma (by purpose I do not write K since this would imply K is dimensionless in usual conventions).

The task is to make gamma equal to whatever you need to make the acceptance probability of the reaction ensemble fit the one of the grand-canonical simulation scheme, see equation (2.60) in my thesis. Next to eq. 2.60 I show how to choose gamma for mimicking a grand-canonical insertion of two monovalent ions (eq. 2.61).

Compare the following example to eq. 2.61: https://github.com/espressomd/espresso/blob/ab38289d308df96c74383d5212b0a4e51b392477/samples/grand_canonical.py#L94 and from the docu http://espressomd.org/html/doc/advanced_methods.html#reaction-ensemble, http://espressomd.org/html/doc/advanced_methods.html#grand-canonical-ensemble-simulation Note, that gamma needs to carry the unit of [sigma^3]^-nubar, where sigma is your unit length of simulation (your simulation box has length x sigma). Note that you need to make sure the unit is correct here! You need to convert whatever concentrations you deal with to 1/sigma^3!

I can see problems coming up easily if you don't follow the given example for the simulation of monovalent salt but want to consider insertions like:

empty <-> 3A +1B

In this case nubar is 4 (and not two like for empty <-> A+B) and then gamma would need to be computed like gamma=c_A_bulk_in_1_div_sigma_cubed^3c_B_bulk_in_1_div_sigma_cubed np.exp(excess_chemical_potential_insertion_of_3A_and_1B(c_A_bulk_in_1_div_sigma_cubed, c_B_bulk_in_1_div_sigma_cubed) / temperature)

Note: excess_chemical_potential_insertion_of_3A_and_1B is a function of c_A_bulk_in_1_div_sigma_cubed, c_B_bulk_in_1_div_sigma_cubed and you have to compute excess_chemical_potential_insertion_of_3A_and_1B for each combination of c_A_bulk_in_1_div_sigma_cubed and c_B_bulk_in_1_div_sigma_cubed you are interested in (for a complex interacting system, there is no way this can be avoided, except in the case of monovalent A and B. See [https://elib.uni-stuttgarI ft.de/bitstream/11682/10848/1/thesis_pdfa.pdf](my phd thesis) for more explainations, epecially pages 21-30, 53-63 for some more background.

If you add ohter ion species C,D;E,... to the reservoir then the excess_chemical_potential_insertion_of_3A_and_1B would not only depend on the concentrations A and B but also C,D,E.... and their excess chemical potentials would also depend on all concentrations.

Also the excess chemical potential depends on the temperature itself, because particles interact differently at different temperatures (for the interaction potentials in question).

I hope this helps you!

Kind regards and a nice weekend

Jonas

dpo854 commented 3 years ago

Thanks for the detailed explanation Jonas. I indeed have gamma in units of 1/sigma^3. And I have determined the chemical potential of both pairs at the desired salt concentration with H+ and OH- concentrations set by the pH of the solution. It looks like I need to do more equilibration runs probably just to reach an initial concentration of ions close to the final density profile.

jonaslandsgesell commented 3 years ago

gamma in units of 1/sigma^3

your gamma has units (1/sigma^3)^nubar, with nubar>1 right? And you test with a system which is filled initally by ions which can all be exchanged with the reservoir (and does not contain particles which cannot be exchanged with the reservoir)? Otherwise you will have the Donnan-Effect and your bulk concentrations won't be reached in the system.

Could you post the values of the excess chemical potentials you determined via widom insertion (at what reservoir concentrations and particle numbers)? Maybe I can see whether the values look sensible.

PS: Probably this is not your main problem, but if you set pH (a measure for the chemical potential of H+ -including the ideal part and the excess part - then the H+ concentration depends on the excess chemical potential of H+ and therefore on all other interactions and concentrations in the system. So it is not quite straight forward to set a pH and salt concentration in a system, chapter 5 in my thesis is about that)

dpo854 commented 3 years ago

In some cases, nubar can be equal to 1. For the exchange of monovalent ions with the reservoir, the chemical potential is cplus cminus exp(excess_chemical_potential/temperature). The same goes for the H+ and OH- pair. Of course, the concentration units are 1/sigma^3. Now, for a salt concentration of 0.01 M and a pH of 4, the chemical potentials for the salt pair and the (H+-OH-) are -0.5484 and -0.5470, respectively. As an additional piece of info, I have a charged wall with dissociable particles in my system. The size of ions is 0.355 nm, which basically is the value of sigma for this case.

jonaslandsgesell commented 3 years ago

Sure, the exchange 1 A2- <-> 2B- would have nubar=1, I intended to exlcude exchanges of single ions from our discussion ;) Also a grand canonical move with only one ion would have nubar=1. However, in this case you deal with a system which is not electro-neutral anymore (because you added 1 ion). Simulations like these typically require quite some care in order not to produce unphysical results.

In the case of the monovalent salt example nubar=2, because empty <-> Na+ + Cl- (and an ion pair is needed for electroneutrality - non-electroneutral macroscopic ssystems are strange).

As an additional piece of info, I have a charged wall with dissociable particles in my system. The size of ions is 0.355 nm, which basically is the value of sigma for this case.

Well in this case i am really, really glad thatyou do not reobtain the (bulk) concentrations you inserted to the grand canonical method because you should not reobtain them due to the influence of the Donnan potential and the walls (which clearly make your system not a bulk system) 🎈 🥳

In order to check whether your grand canonical setup has no errors you need to simulate a bulk system (no walls, only exchangeable ions, no other particles like counterions of the wall etc.) and see whether you obtain the concentrations you set. This has to be true (with small deviations due to sampling). If it is true, you checked your grand canonical part of the simulation successfully and you can trust it.

Other stuff which could go wrong is the electrostatic integration/energy computation in the slab etc. This you would want to check separately with the method of your choice.

PS: > -0.5484 and -0.5470 [kT] from a first glance this looks reasonable for the excess chemical potenital of a system of only monovalent ions at kT=1 and bjerrum length 0.71nm (see figure 2.3 in my thesis). For multivalent ions, I suggest to make the check above!

PPS: the difference between the chemical potential and the excess chemical potential is the ideal/kinetic part kTln(c:i*lambda_B,i^3) which is typically the dominating part. That is why in chemistry it is quite common to formulate the law of mass action with the concentrations (neglecting excess contributions in the activity coefficients)

dpo854 commented 3 years ago

That's a good benchmarking suggestion. Will surely put my GC steps through this test.

kosovan commented 3 years ago

Regarding brute-force sampling a rare event, please note the following:

I agree with Jonas that rare insertions at very close distances spoil the averaging. If you look at the time series, you can notice spikes in it which correspond to these insertions.
The short-range potentials used in CG models do not provide realistic description at short distances r<<sigma, so you are probably not interested in accurately calculating the short-range contribution at these distances.
When you combine short-ranged repulsion with Coulomb attraction at $r\to 0$, you get $\infty - \infty$ which can yield any result if you consider rounding errors.

Few years ago, with Tobi Richter we tested the following workaround:

If you use LJ potential in MD simulation, your particles never come closer than $r\approx 0.9\sigma$. Therefore, for short distances you can safely define $U(r<0.8\sigma)=+\infty$ without affecting the result.
Use Widom to explicitly sample distances $r>0.8\sigma$ where the short-range interaction is relevant to compute the excess chemical potential. This does not suffer from the problems with too-close insertions.
To this contribution you should add excess chemical potential of hard spheres with radius $R=0.8\sigma$ computed analytically from your favourite equation of state of hard spheres, e.g. Carnahan-Starling.

This approach gave us results which were very consistent with the chemical potential calculated from the grand-canonical insertion of ion pairs even at rather high densities. I am sorry, I do not remember the exact numbers because we never published these results. But I encourage you to try it out.

jonaslandsgesell commented 3 years ago

When you combine short-ranged repulsion with Coulomb attraction at $r\to 0$, you get $\infty - \infty$ which can yield any result if you consider rounding errors.

I never observed arbitrary results in Espresso but the interaction energy was always reported correctly as diverging to plus infinity (within numerical precision). Do you have a minimal working example for this case? I would consider this a bug if it really happens and the r^-12 term was not dominating over the r^-1 coulomb term in taking the limit r to 0 https://www.wolframalpha.com/input/?i=Lim+r+to+0+%28A%2Fr%5E12-B%2Fr%5E6-C%2Fr%29

kosovan commented 3 years ago

Do you have a minimal working example for this case? I don't have a minimal working example. Probably you are right hat Espresso will always give you +infty. However, this does not change anything on the rest of my comment.

The term $r^12$ in LJ potential is not based on physical significance - it was introduced for mathematical convenience. Therefore, it makes no sense to use sophisticated methods and accurately sample its contribution. Any other diverging potential will do the job of representing steric repulsions at $r \to 0$, so you can arbitrarily choose a more convenient one. If you do it only in the region which is inaccessible to thermal fluctuations, you will not notice a difference.

jonaslandsgesell commented 3 years ago

If you do it only in the region which is inaccessible to thermal fluctuations, you will not notice a difference.

I completely agree! One can do the proposed exchange of the short ranged diverging potential with neglegible error in this case.

However, as far as I recall (from long time ago), the jumps in the chemical potential were because of insertions around 0.8 sigma <r < 4sigma (?) So it was not the short range diverging potential but the attractive part of the coulomb interaction which made the "jumps" in the excess chemical potential estimate happen. Can we check that?

PS: if this is true, then the repulsive "soft" core potential is not the probelm and jumpiness of the excess chemical potential estimate cannot be removed by this hard core interaction trick.

jonaslandsgesell commented 3 years ago

I made a quick run which supports my observation. In the log file, before each reported pair (value of measured mu_ex, err_mu_ex) there is the closest distance to each inserted test particle reported (for each test particle separately). In the log the first jump occurs after inserting an ion with distance 1.37116 and another ion with distance 5.31111 to an existing particle: ... (-0.8346797576640681, 0.3069210010050128) d_min 1.37116 d_min 5.31111 (-2.298778979964505, 0.7716876636582033)

another jump occurs here:

(-1.28650655997571, 0.09712397377552845) d_min 1.72205 d_min 8.404 (-1.3244043386740338, 0.100548296368586)

another one here

(-1.2820831487646223, 0.07491243174232182) d_min 2.28557 d_min 1.35472 (-1.3372468339670849, 0.08888452185707087)

I have the impression that insertions where particle do overlap however are not the reason for jumps in the excess chemical potential estimate:

The first occurance is here: (-1.6729595512991167, 0.28767162934202123) d_min 0.681229 d_min 6.04046 (-1.6665286609688261, 0.28773797994322386) which barely changes the excess chemical potential estimate

log_monovalent_ions_cbulk_0.0005_in_1_div_sigma_cubed_bjerrum_length_6_sigma_kT_1.txt

Therefore, I do not expect that making the soft core LJ potential behave like a hard core interaction for distances smaller than 0.8 sigma would make the estimate of the excess chemical potential less susceptible to "jumps".

jonaslandsgesell commented 3 years ago

I think this issue can be closed. The vanilla Widom insertion estimate of the excess chemical potential is to my understanding inherently sucseptible to jumps due to "rare" events. Widom insertion needs to be run long enough for stable results.

dpo854 commented 3 years ago

For sure the classical widom insertion method is accurate enough for what you want to do. In order to be on the safe side, just run a grand canonical simulation with the excess chemical potential for a given conce traction you determined as input and some random initial number of particles in the box. Then check whether the concentrations you want to achieve are actually achieved in the simulation box after equilibration.

Unfortunately, the number of particles in the system comes out to be a bit greater than the number corresponding to the desired concentration. To be more specific, I did the simulation with a volume corresponding to 1000 particles for a given concentration and the average number of particles in the system is around 1175 with the imposed chemical potential. I believe a finite value of the exclusion radius might have played a role in that since a larger exclusion radius would equate to more acceptances.

jonaslandsgesell commented 3 years ago

Unfortunately, the number of particles in the system comes out to be a bit greater than the number corresponding to the desired concentration. To be more specific, I did the simulation with a volume corresponding to 1000 particles for a given concentration and the average number of particles in the system is around 1175 with the imposed chemical potential. I believe a finite value of the exclusion radius might have played a role in that since a larger exclusion radius would equate to more acceptances.

I think you refer to the "empty" box test? ("empty" here means that you setup the box with some number of exchangeable ions (and nothing else) for optimal tuning of P3M and then you see whether the average concentration you obtain in the end is the concentration you set in the grand canonical scheme) The concentration you encounter is 17.5% increased compared to what you expect (1175/1000). What is the average concentration fluctuation? How long did you equilibrate your system? Can you check equilibration e.g. make sure the particle numbers do not have a trend anymore. Maybe it is worth to check the grand canonical scheme with a box which has a smaller volume, there particle numbers will be smaller and the computations can be performed faster -- allowing for longer equilibration times. I ended up doing on the order of 1.5 million MC moves for reactions/particle insertions for some simulations prior to starting a run in order to ensure equilibration and be on the safe side. The exclusion radius itself won't affect the equilibrium state (because detailed balance is obeyed) but an exclusion radius too big can make your acceptance ratio very small. This means you should try to choose an exclusion radius as small as possible (which your MD integrator still can handle, typical values which are good are 0.9 sigma).

kosovan commented 3 years ago

Unfortunately, the number of particles in the system comes out to be a bit greater than the number corresponding to the desired concentration. To be more specific, I did the simulation with a volume corresponding to 1000 particles for a given concentration and the average number of particles in the system is around 1175 with the imposed chemical potential. I believe a finite value of the exclusion radius might have played a role in that since a larger exclusion radius would equate to more acceptances.

This is the expected correct result because the excess chemical potential is negative, Fig. S3 of the supporting information for this paper

The full chemical potential is the sum of reference, ideal and excess terms, mu = mu_ref + mu_ideal + mu_ex If mu_ex < 0, then mu_ideal = kT ln (c / c_ref) = mu - mu_ex -mu_ref > mu - mu_ref

Exclusion radius plays no role - it only affects the efficiency (how quickly you converge to the result) but not the value to which you converge.

kosovan commented 3 years ago

The system has two pairs: the first, cation=+3 and anion=-1, and the second, H+ and OH-. The jump is highlighted in bold.

@dpo854 Perhaps the jump might be related to the fact that you happened to insert the ion close to a 3+ ion that you have in your system. This might be indeed a rare event if you have a low number of 3+ ions. Could you please post your simulation script, so that we could look at the parameters and the whole setup?

jonaslandsgesell commented 3 years ago

This is the expected correct result because the excess chemical potential is negative

It depends on how the test is performed right? If you set certain bulk concentrations and the corresponding excess chemical potential for this bulk concentration in the grand canonical insertion method (i.e. you set a chemical potential), then the expected behavior is to obtain the bulk concentrations in the "empty box" test as average concentrations in the box.

dpo854 commented 3 years ago

Unfortunately, the number of particles in the system comes out to be a bit greater than the number corresponding to the desired concentration. To be more specific, I did the simulation with a volume corresponding to 1000 particles for a given concentration and the average number of particles in the system is around 1175 with the imposed chemical potential. I believe a finite value of the exclusion radius might have played a role in that since a larger exclusion radius would equate to more acceptances.

I think you refer to the "empty" box test? ("empty" here means that you setup the box with some number of exchangeable ions (and nothing else) for optimal tuning of P3M and then you see whether the average concentration you obtain in the end is the concentration you set in the grand canonical scheme) The concentration you encounter is 17.5% increased compared to what you expect (1175/1000). What is the average concentration fluctuation? How long did you equilibrate your system? Can you check equilibration e.g. make sure the particle numbers do not have a trend anymore. Maybe it is worth to check the grand canonical scheme with a box which has a smaller volume, there particle numbers will be smaller and the computations can be performed faster -- allowing for longer equilibration times. I ended up doing on the order of 1.5 million MC moves for reactions/particle insertions for some simulations prior to starting a run in order to ensure equilibration and be on the safe side. The exclusion radius itself won't affect the equilibrium state (because detailed balance is obeyed) but an exclusion radius too big can make your acceptance ratio very small. This means you should try to choose an exclusion radius as small as possible (which your MD integrator still can handle, typical values which are good are 0.9 sigma).

Yes, this is what I am referring to. Let me put your suggestion to the test and see if it can improve the agreement.

dpo854 commented 3 years ago

I did a test with 100 particles in the box with 1 million MC moves before taking the data. The average number of particles and the standard deviation are 118.167 and 0.07167, respectively. The statistics were collected over 10000 cycles where each cycle consists of 20 reactions and 200 MD runs.

jonaslandsgesell commented 3 years ago

I think this deviation is too big.

Can you please post your simulation script (grand canonical and Widom insertion, both) and the parameters with which you called them?

In the mean time I suggest to run the same test with monovalent ions. There we proved that this test converges.

Thank you!

dpo854 commented 3 years ago

This is Widom.

import numpy as np
import argparse
import statistics 

import espressomd
from espressomd import reaction_ensemble
from espressomd import electrostatics
from espressomd.io.writer import vtf
import datetime

required_features = ["P3M","LENNARD_JONES"]
espressomd.assert_features(required_features)

sigma=0.355
cs_bulk=0.01
cs_bulk = (cs_bulk)*0.6023*(sigma**3)
N0 = 100
qcat=1.0
qana=-1.0
box_l = (N0 / cs_bulk)**(1.0 / 3.0)
print("boxl ",box_l)
volume=(box_l)**3
pHbulk=4.0
pH=pHbulk
pOH=14-pH
cHplus=10**(-pH)
cOHminus=10**(-pOH)
noHplus=int(cHplus*0.6023*(sigma**3)*volume)
noOHminus=int(cOHminus*0.6023*(sigma**3)*volume)
nopositive=int(N0)
nonegative=int(N0)*(round(qcat))

if(pH<7):
  surplus=noHplus-noOHminus
  extraana=surplus/(-round(qana))
  remainder=surplus-extraana*(-round(qana))
  nonegative=nonegative+round(extraana)
  noOHminus=noOHminus+round(remainder)
  surplus=(noHplus+nopositive*round(qcat)-noOHminus+nonegative*round(qana))
  print("surplus ",surplus)
  if(surplus>0):
     noOHminus=noOHminus+surplus
  else:
     noHplus=noHplus+abs(surplus)
else:
  surplus=noOHminus-noHplus
  remainder=surplus%(round(qcat))
  extracat=surplus/(round(qcat))
  nopositive=nopositive+round(extracat)
  noHplus=noHplus+round(remainder)
  surplus=(noHplus+nopositive*round(qcat)-noOHminus+nonegative*round(qana))
  print("surplus ",surplus)
  if(surplus>0):
     noOHminus=noOHminus+surplus
  else:
     noHplus=noHplus+abs(surplus)

system = espressomd.System(box_l=[box_l, box_l, box_l])
np.random.seed(seed=42)
system.time_step = 0.01
system.cell_system.skin = 0.4
temperature = 2.4935

# type 0 = MX
# type 1 = Mx+
# type 2 = Az-
# type 3 = HOH
# type 4 = H+
# type 5 = OH-

for i in range(nopositive):
    system.part.add(id=i, pos=np.random.random(3) * system.box_l, type=1, q=qcat)
for i in range(nopositive, nopositive+nonegative):
    system.part.add(id=i, pos=np.random.random(3) * system.box_l, type=2, q=qana)
for i in range(nopositive+nonegative, nopositive+nonegative+noHplus):
    system.part.add(id=i, pos=np.random.random(3) * system.box_l, type=4, q=+1.0)
for i in range(nopositive+nonegative+noHplus,nopositive+nonegative+noHplus+noOHminus):
    system.part.add(id=i, pos=np.random.random(3) * system.box_l, type=5, q=-1.0)# all particles added

lj_eps = 1.0
lj_sig = 1.0
WCA_cut=2.0**(1.0/6.0)
lj_cut=WCA_cut*lj_sig
types = [0, 1, 2, 3, 4, 5] # 0 MX, 1 Mx+, 2 X-, 3 HOH, 4 H+, 5 OH-
for type_1 in types:
    for type_2 in types:
        system.non_bonded_inter[type_1, type_2].lennard_jones.set_params(epsilon=lj_eps,sigma=lj_sig, cutoff=lj_cut, shift=0)

p3m = electrostatics.P3M(prefactor=5.081, accuracy=1e-3)
system.actors.add(p3m)
p3m_params = p3m.get_params()
for key, value in p3m_params.items():
    print("{} = {}".format(key, value))

warm_steps = 50
warm_n_times = 50
min_dist=0.9*lj_sig

system.integrator.set_steepest_descent(f_max=0, gamma=1e-3,
                                       max_displacement=0.01)
i = 0
while system.analysis.min_dist() < min_dist and i < warm_n_times:
    print("minimization: {:+.2e}".format(system.analysis.energy()["total"]))
    system.integrator.run(warm_steps)
    i += 1

print("minimization final: {:+.2e}".format(system.analysis.energy()["total"]))
print()
system.integrator.set_vv()

system.thermostat.set_langevin(kT=temperature, gamma=1.0, seed=42)

widom = reaction_ensemble.WidomInsertion(
    temperature=temperature, seed=77)

insertion_reaction_id1 = 0
insertion_reaction_id2 = 1
widom.add_reaction(reactant_types=[],
                   reactant_coefficients=[], product_types=[1, 2],
                   product_coefficients=[1, round(qcat)], default_charges={1: qcat, 2: qana}) # for the reaction MXn -> Mn+ + n X- 
widom.add_reaction(reactant_types=[],
                   reactant_coefficients=[], product_types=[4, 5],
                   product_coefficients=[1, 1], default_charges={4: +1.0, 5: -1.0}) # for the reaction HOH -> H+ + OH- 

print(widom.get_status())
system.setup_type_map([0, 1, 2,3,4,5])
fp = open('trajectory.vtf', mode='w+t')
vtf.writevsf(system, fp)
vtf.writevcf(system, fp)
clock_time1=datetime.datetime.now()
fin1=open("potential.txt","w")
n_iterations = 1000
pressurelist=[]
system.integrator.run(steps=100000)
for i in range(n_iterations):
    for j in range(100):
        widom.measure_excess_chemical_potential(insertion_reaction_id1)
        widom.measure_excess_chemical_potential(insertion_reaction_id2)
    system.integrator.run(steps=1000)
    pressure=system.analysis.pressure()
    press=pressure["total"]
    pressurelist.append(press)
    vtf.writevcf(system, fp)
    if i % 20 == 0:
        print("mu_ex_pair 1 ({:.4f}, +/- {:.4f})".format(
            *widom.measure_excess_chemical_potential(insertion_reaction_id1)))

        print("mu_ex_pair 2 ({:.4f}, +/- {:.4f})".format(
            *widom.measure_excess_chemical_potential(insertion_reaction_id2)))

print("excess chemical potential for an ion pair MXn ",
      widom.measure_excess_chemical_potential(insertion_reaction_id1))

print("excess chemical potential for an ion pair HOH ",
      widom.measure_excess_chemical_potential(insertion_reaction_id2))
clock_time2=datetime.datetime.now()
print("Time elapsed ",clock_time2-clock_time1)
meanpress=statistics.mean(pressurelist) 
stdevpress=statistics.stdev(pressurelist) 
fp.close()
print("mean press ",meanpress, "stdev press", stdevpress)

The salt concentration is 0.01 M and the pH, 4. The ions are monovalent. So, there is greater certainty about the chemical potential computation being correct since it's close to Debye-Huckel predictions, which works well for monovalent ions.

dpo854 commented 3 years ago

The lines are getting distorted for some reason. I am sending the codes over by email.

jonaslandsgesell commented 3 years ago

I spotted a unit conversion problem in this line: cs_bulk=(args.cs_bulk)**0.6023*(sigma**3) # converting it into 1/sigma**3 which should rather read: cs_bulk=(args.cs_bulk)*0.6023*(sigma**3) # converting it into 1/sigma**3 This was the main problem.

Other hints from our private communication are:

RE.add_reaction(gamma=cHplus*cOHminus * np.exp(excess_chemical_potential_pair_HOH/temperature),reactant_types=[], reactant_coefficients=[], product_types=[types["H+"], types["OH-"]],product_coefficients=[1, 1], default_charges={types["H+"]: +1, types["OH-"]: -1})

Please note that like this H+ and OH- can only fluctuate together (e.g. go 1 up in particle number together or go 1 down in particle number together, because you describe insertion/deletion of pairs). If you first have to find the correct amount of H+/OH- in your box, you need to have them fluctuate independently (like I did in my phd thesis and the paper https://pubs.acs.org/doi/full/10.1021/acs.macromol.0c00260). You can do this if you come up with inserting (H+ + salt anion) and (OH- + salt cation)... like it is outlined in section "5.3.2 Particle Exchanges" in my phd thesis https://elib.uni-stuttgart.de/bitstream/11682/10848/1/thesis_pdfa.pdf

Please run the Widom insertion much longer than it is done in the script right now. For me the excess chemical potential still changed in the second digit after the point. The longer, the better. Also the Debye-Hueckel approximation will break down for high electro static couplings or multivalent ions...
You might also consider using pint for making unit conversions: see e.g. here https://github.com/espressomd/espresso/blob/python/doc/tutorials/constant_pH/constant_pH.ipynb and https://pint.readthedocs.io/en/stable/

PS: if all your ion pairs have the same interaction potentials and parameters you only have to determine one excess chemical potential for one ion pair and reuse it.

dpo854 commented 3 years ago

Thanks a lot Jonas. I could never imagine that a typo would be at the heart of all these problems. I will try to be more careful not to make such errors in the future. And also a lot of thanks for the resources for going deep into the inticracies of pH-dependent simulations.

jonaslandsgesell commented 3 years ago

Hey Debadutta,

You are welcome! I am glad I could help you. I am still graceful for how other people also helped me :)

dpo854 commented 3 years ago

Just an update on the GC simulation. Now, the average number of particles in the system comes out to be 99.276 with a standard deviation of 0.04. I guess that sorts out the problem.

espressomd / espresso

Variance in chemical potential computation #4003