thank you for your help. I am writing to ask whether I am missing something, or the gmx bar function is not printing any DG contribute from the mutation of masses: I am trying one “real” case (see first column of every .xvg file in the get_results_proper folder), and one dummy case (first column scores -1000 everywhere), but still having same result: 0.00 kJ/mol.
As a second topic, I would like to now if the case of couple-intramol = yes keyword actually does not count at all the 1-2 and 1-3 charge and vdWs interactions within the molecule (i.e. they are not contributes in dH/dlambda), as they’re not even applied given the bonded parameters parametrization.
I haven’t perturbed masses myself, but I don’t think it is accounted for when calculating dH/dL. I don’t think it should contribute either, but I may be wrong.
couple-intramol = yes means that the intramolecular interactions of the decoupled molecules are also being decoupled (except 1-4 interactions). It does not add any interactions that were not there in a non-alchemical FEP simulation.
Mass changes should contribute to dH/dlambda. There is code that adds these contributions, also to delta H for foreign lambda’s. I assume this should work, but there can always be some bug in the logic.
I am sorry I am late on replying, and moreover partial in my question. I really meant “add” as “numerically counted” in the calculation of the free energy, I read my second paragraph now and notice it is unclear and might be read as literally “adding interactions” in the MD run.
I took some time to better assess some details, I am attaching them below since you always helped me a lot, hope at least this provides some clarifications/curiosities in turn; thank you.
Kindly thank you both. I wanted to carry out some additional verifications, thus I am providing 4 different examples. I agree with the observation that the kinetic part should be included due to the kinetics term dependence on m(lambda), but unfortunately I had to correct my reply of this morning and tell that I don’t think this counted. The following will tell if I am wrong; I am writing something you can see in the README in the link below, the names reported refer to the corresponding directory:
‘H2O_to_NA’ shows that couple-intramol=yes does not consider up to the 1-3 intramolecular q and vdW interactions (it is actually an NA —> H2O, but not important);
‘q_14_counted’ shows that couple-intramol=yes considers 1-4 intramolecular interactions, scales them and count the contribute in DeltaG (DeltaA);
‘morphe_mass’ shows that mass mutation is the detected, but not counted in the result (bug? If not, Prof. Hess is right, I am missing the reason why);
‘q_and_mass’ shows that, if in the topology the B state is with different charges, but the coul-lambdas = array_of_zeros , the code correctly avoids computing the q contribute to DG, but detects of the q change;
The transformations are immediately understood by looking at the corresponding topol.top (merely at the [ atoms ] directive) and md.mdp files (merely at the lambda series for q, vdw, masses).
I am sorry, I may sound naif, but actually I want to point out that what worries me is that the ‘morphe_mass/get_results’ directory shows that the xvg files contained manifest large dh/dl values related to masses, but then the ‘results.txt’ file shows just a sequence of zeros for DeltaGs. At the same time, those xvg files do not show any influence of mass variation on the “s5” and “s7” columns within themselves, even though the mdp files adopted showing that mass_lambdas is set as a non-zeros sequence.
Sorry I am also for the 1-3 related topic question, which I was aware of, but since this putable problem I think I have with masses, I wanted to carry out further sanity checks when plugging free energy methodology, as my project depends on that.
Thank you for your time, maybe the ‘morphe_mass’ example is just an incorrect interpretation of a beginner?
Ah, I see the issue now. I was looking at dH/dlambda, as in the title of your post, that is present. But the foreign lambda terms are all zero and that is what gmx bar uses.
thanks to all of you, instead. I am sharing another bug about the free energy code, regarding the decoupling from the environment of a charged solute (i.e., switching off its charges and vdWs): Free energy, annihilation: charged ligand - #4 by Jacopo
In that very case, the problem concerns the fact that no warning about net charge of the box in any intermediate (and the ending) state is ever raised. I shared examples in there and to be crystal clear and not waste your time.
Hoping this helps, I am looking forward to any news, wish you a good day.