Hello,
I have a standard periodic box protein ligand system containing 9 Na+ counterions for neutralisation purposes.
The ligand contains 2 phosphate groups which interact with the protein and throughout MD simulations several Na+ ions can come and stick to these phosphate groups which is not ideal/wanted.
I've therefore been working on trying to place harmonic restraints to prevent the sodium ions from coming too close to the phosphate groups.
I've tried a few different schemes, in each case, each ion has its own restraint and it is to a group of backbone atoms whose combined center of mass is close to both phosphate groups. I've also tried more direct restraints between the phosphate groups and each ion (I wanted to avoid restraints on the phosphate groups to avoid disturbing it, but that's off the topic for this email).
I see in the log files the restraints are read in and the restraints themselves are working/doing what is expected most of the time. The issue is after some time (in the 10s of ns), I can observe a sodium ion interacting with a phosphate group after post-processing the trajectory. When I look at the DISANG file (containing the calculated restraint distances, the distances are always above the thresholds). Sometimes, shortly after this the simulation crashes (presumably due to the high forces generated by the restraints after the system images the ion, resulting in the ion being very close to the ligand).
My thinking/understanding is this is an imaging/periodicity issue in which the restraint distances are not correctly imaged at every frame meaning the ion can in some cases come very close to the ligand without the restraint activating. Does that sound right?
Looking through the forum I seen suggestions like setting "nscm" to be greater than the timestep. I tried that but it unfortunately did not work for me. I've also played with "iwrap".
I also tried to use plumed to define the restraints but ultimately had the same issue.
Does anyone have any suggestions on what could be done/tried? My only thought at the moment is to just remove the counter ions and use the neutralising plasma approach.
Thanks for reading my question.
Kind regards,
Rory
Below is an example production md input file and restraint file, happy to provide more inputs/info if useful of course
############# production file ###########
400 ns NPT MD production
&cntrl
imin=0, irest=1, ntx=5, iwrap=0,
nstlim=200000000, dt=0.002,
ntpr=5000, ntwx=5000, ioutfm=1, ntwr=50000,
ntc=2, ntf=2, cut=8.0,
ntb=2, ntp=1, pres0=1.0, taup=1.0,
ntt=3, tempi=298.0, temp0=298.0, gamma_ln=5.0, ig=-1,
nmropt=1, nscm=200000001,
/
&wt type='DUMPFREQ', istep1=5000,
/
&wt type='END'
/
DISANG=distrestraint.RST
DUMPAVE=dist_vs_t.txt
############# restraint file ###########
&rst
r1=18,r2=20,r3=45,r4=47,
rk2=10.0,rk3=0.0,
iresid=0,
iat=-1,4003
igr1=173,175,183,184,756,758,774,775,2002,2004,2010,2011,2012,2014,2032,2033,2650,2652,2662,2663
/
&rst
r1=18,r2=20,r3=45,r4=47,
rk2=10.0,rk3=0.0,
iresid=0,
iat=-1,4004
igr1=173,175,183,184,756,758,774,775,2002,2004,2010,2011,2012,2014,2032,2033,2650,2652,2662,2663
/
....
[repeated for each ion numbered between 4003 and 4014]
...
&rst
r1=18,r2=20,r3=45,r4=47,
rk2=10.0,rk3=0.0,
iresid=0,
iat=-1,4014
igr1=173,175,183,184,756,758,774,775,2002,2004,2010,2011,2012,2014,2032,2033,2650,2652,2662,2663
/
--
Rory Crean
När du har kontakt med oss på Uppsala universitet med e-post så innebär det att vi behandlar dina personuppgifter. För att läsa mer om hur vi gör det kan du läsa här: http://www.uu.se/om-uu/dataskydd-personuppgifter/
E-mailing Uppsala University means that we will process your personal data. For more information on how this is performed, please read here: http://www.uu.se/en/about-uu/data-protection-policy
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Apr 25 2023 - 03:30:02 PDT
