normal Running alchemical amino acid transformation with Martini 3

  • wuzhiyi
  • wuzhiyi's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
3 years 7 months ago #8734 by wuzhiyi
I’m interested in running amino acid alchemical transformation with martini3. For example to change isoleucine to valine,
I change the atom
467 SC1 189 ILE SC1 467 0.0000 72.0 SC2 0.0 72.0; H
and the bond contraint
466 467 1 0.34100 0.29200 ; ILE > VAL
However, I noticed that when I try to run the simulation with bonded-lambdas, I got the error

Fatal error:
Function type Restricted Angles does not support currentely free energy
calculations

This seems to related to the type 10 angle used which is not perturbed and it seems that having the type 10 angle alone will trigger this fatal error and one cannot run any alchemical transformation with change in bonded terms.

Please Log in or Create an account to join the conversation.

More
3 years 7 months ago #8735 by riccardo
True, Gromacs doesn't support free energy (FE) calculations with a few function types including the Restricted Bending (ReB) potentials.

The workaround is to substitute that potential with a classic harmonic or cosine harmonic potential (function types 1 or 2). They will be similar enough close to the equilibrium position so as to not to really make a difference to your FE calcs - see also the plots in the paper which introduces the ReBs pubs.acs.org/doi/abs/10.1021/ct400219n .

The ReBs are there for numerical stability purposes - see: pubs.acs.org/doi/abs/10.1021/ct400219n . Even without them, you should still be able to run the FE calcs; if your system is big (i.e., contained many of such ReBs), it might be that you experience a few simulation crashes, but just re-starting the simulations should do.

Hope this helps!

Please Log in or Create an account to join the conversation.

  • wuzhiyi
  • wuzhiyi's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
3 years 7 months ago #8736 by wuzhiyi
Hi riccardo,
Thank you for your reply. I have tried to change the angle type from 10 to 2. However, I noticed that I got LINC error due to probably too low the restraint strength of the SBB and BBS angle.
I'm thinking of increasing the restraint constant to a value higher than the current 5. The martini2 uses a restraint constant of 25 paired with angle type 2 but has only one pair of restraint for each residue. I wonder if you have any suggestions with regard to how high should these restraint constants be?

Please Log in or Create an account to join the conversation.

  • wuzhiyi
  • wuzhiyi's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
3 years 7 months ago #8737 by wuzhiyi
Hi riccardo,
Thank you for your reply. I have tried to substitute the SBB and BBS angle from 10 100 5 to angle type 1 or 2 and using a restraint strength of up to 25. However, it doesn't seem that the angle restraint allows a stable simulation where I always got LINC error.

Please Log in or Create an account to join the conversation.

More
3 years 7 months ago #8739 by riccardo
Hi wuzhiyi and sorry for the slow follow-up. How many of such ReBs were present in the itp? I was hoping it would run and just crash a few times. Is the system properly equilibrated, i.e., does it run nicely if the ReBs are in place?
Regarding the question about the force constants, I'm not sure, you could definitely try to play around with the (i.e., increase them) to see whether they help or not with the numerical instabilities.
One last option would be to lower your time step just for the free energy calculations.

Please Log in or Create an account to join the conversation.

  • wuzhiyi
  • wuzhiyi's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
3 years 7 months ago #8740 by wuzhiyi
Hi riccardo,
Thank you for the help.
I'm simulating a protein with ~200 residues so around 300~400 ReBs
with the format of

2 1 3 10 100 5 ; M(1)-M(1)-N(2) SBB
1 3 4 10 100 5 ; M(1)-N(2)-N(2) BBS

If the ReBs are there the system is very stable. It seems that increasing the force constant will help. However, I have increased it to 25 and it still gives linc error after running for a while. I'm not sure if increasing it to a higher value is ok.
I'm currently using a timestep of 20fs.

Please Log in or Create an account to join the conversation.

More
3 years 7 months ago #8741 by riccardo
Alright, that's quite a number of ReBs.
Well, what the implications of using a higher force constant for these angles are depends on what you are interested in doing with these calculations. You should try to assess the impact on the quantity you want to compute of having different force constants on those angles.

Now, one option is that you write a little script that resubmits the simulation whenever it crashes (using the Gromacs checkpoint file). Of course, this makes sense if your simulation runs for a "decent" amount of time.

Finally, as hinted at before, you could choose to lower the time step. This should improve numerical stability. In this way, you can perhaps have *numerically stable but slower* simulations. Perhaps using 15 fs or 10 fs already works for you.

Please Log in or Create an account to join the conversation.

  • wuzhiyi
  • wuzhiyi's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
3 years 7 months ago #8742 by wuzhiyi
Hi riccardo,
Thank you for your help.
I wonder do you know when will martini3 available? If it will take a year or two, I will try to see which of the methods can solve my question. otherwise, I will probably just use non-alchemical methods.

Please Log in or Create an account to join the conversation.

Time to create page: 0.103 seconds