Inquiry about estimated errors.

[bowensha-tudelft]

Hi. When I was using pAPRika, I am always confused by the error it gives. What can I do to reduce the error? I have tried to increase number of windows and elongate time scale of my simulations. But it seems they do not work very much.

[jaketanderson]

Hi Bowen, what's the particular error you see?

[bowensha-tudelft]

The error given by paprika. Maybe it is called SEM? I tried to read the original paper but could not understand it well.

[jset231]

@bowensha-tudelft, Can you print or take a snapshot of the error you're getting? Maybe that would help us debug the issue

[slochower]

I think @bowensha-tudelft may be asking about reducing the uncertainty in the estimate of dG. If that's the case, running longer in each window, using more windows, tweaking the force constants are all options. In general, the uncertainty in the estimate will go down with more sampling. There may be other sources of error, however, like the force field used to model the host and guest molecules.

[bowensha-tudelft]

@bowensha-tudelft, Can you print or take a snapshot of the error you're getting? Maybe that would help us debug the issue

I am trying to run notebook-01 in paprika tutorials. I changed the host and guest to beta-cyclodextrin and hexanoate. I followed settings in the paper (https://pubs.acs.org/doi/10.1021/acs.jctc.5b00405). But the result is quite different.

I tried to make a graph of energy during the process (below). I use command like "plt.plot(attach_fractions,free_energy.results["attach"]["ti-block"]["fe_matrix"][0,:])". I am not sure if i use the right energy.

When I in crease the force constants of dihedrals within beta-cyclodextrin, the energy will rise in attach phase (blue lines). Now the maximum is about 25 and if the force constants are larger, it will be like 75 but it does not fall in release phase. I guess I did something wrong. And sometimes the calculated binding affinity of bcd and hex is less than -100 kcal/mol. That definitely means there are many errors in my settings.

Here is my settings for constranints:

They are 14 dihedrals in beta-cyclodextrin, also mentioned in the paper.

[jset231]

Hmm, that is interesting; it could result from many different things. Just to make sure, are you getting a delta G (bind) around -100kcal/mol? That seems quite large and probably has something to do with the restraints. In explicit solvent, you should get closer to -10 kcal/mol.

Perhaps the atom mask is not specified correctly. As a sanity check, print out the atom indices of the DAT restraints and check if they correspond to the correct atoms.

for res in guest_restraints:
    print(res.index1, res.index2, res.index3, res.index4) 

Also for the target and force constant values, make sure you set the correct units with OpenFF-Units. For example:

from openff.units import unit as off_unit

r.attach["target"] = -112.5 * off_unit.degrees
r.attach["fc_final"] = 6.0 * off_unit.kilocalorie_per_mol / off_unit.degrees**2