prody / coMD

This package runs collective molecular dynamics simulations, where collective motions from anisotropic network model normal mode analysis drive all-atom dynamics through targeted molecular dynamics simulations and energy minimization
http://prody.csb.pitt.edu/comd/
MIT License
1 stars 6 forks source link

coMD - walker 2 TMD problem in comd.tcl #62

Closed abvl closed 5 years ago

abvl commented 5 years ago

I would like to use coMD for pathway analysis but noticed a problem in the comd.tcl files on both the prody website and github pages...

using the comd tutorial files, if I run with either 1ake or 4ake on their own (ie same filename for initial and final) everything is fine.

However, when I use two different files to generate a pathway it fails to run TMD on walker 2 and it seems to be because the code is trying to run namd using the wrong .conf file...

line 1317 of comd.tcl currently writes to 'pro.conf' but further down at line 1377 it tries to run namd using 'min.conf' - changing this to 'pro.conf' by analogy to the sections for Walker 1 allows everything to run and I can see the process go through several cycles fine.

However, I'm not clear about what parameters to use for ANM and TMD to complete a successful run using the tutorial files - it crashes every time after a few cycles. Looking at the '...anmmc_log.txt' files shows that when things crashe the 6th column contains rapidly increasing values that end up as 'inf' .

Could you please post a set of values for the ANM and TMD options that should work with the tutorial ?

Thanks,

Andrew Beavil

jamesmkrieger commented 5 years ago

Hi Andrew,

You can find values in the original paper on coMD with adenylate kinase: https://www.ncbi.nlm.nih.govpubmed/24094405/

Thanks for reporting your problem. I'll have a look into it and see if I can work out what's going wrong.

Best wishes James

On Nov 28, 2019, at 5:31 AM, abvl notifications@github.com wrote:

 I would like to use coMD for pathway analysis but noticed a problem in the comd.tcl files on both the prody website and github pages...

using the comd tutorial files, if I run with either 1ake or 4ake on their own (ie same filename for initial and final) everything is fine.

However, when I use two different files to generate a pathway it fails to run TMD on walker 2 and it seems to be because the code is trying to run namd using the wrong .conf file...

line 1317 of comd.tcl currently writes to 'pro.conf' but further down at line 1377 it tries to run namd using 'min.conf' - changing this to 'pro.conf' by analogy to the sections for Walker 1 allows everything to run and I can see the process go through several cycles fine.

However, I'm not clear about what parameters to use for ANM and TMD to complete a successful run using the tutorial files - it crashes every time after a few cycles. Looking at the '...anmmc_log.txt' files shows that when things crashe the 6th column contains rapidly increasing values that end up as 'inf' .

Could you please post a set of values for the ANM and TMD options that should work with the tutorial ?

Thanks,

Andrew Beavil

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

jamesmkrieger commented 5 years ago

Oh I should also note that anmmc.py has default values. That will be used if you leave the boxes blank.

jamesmkrieger commented 5 years ago

I am now using chain A from 1ake and chain A from 4ake, and running a test run with default values and I'll let you know if it works.

jamesmkrieger commented 5 years ago

It worked for me with 6 cycles after making the fix you found. Thanks for reporting it. I will now merge it so other users can get that too.

Best wishes James

jamesmkrieger commented 5 years ago

What values did you use that caused problems? Could it be that you forgot to select only chain A?

abvl commented 5 years ago

Dear James,

thanks for responding so quickly...

on reflection, I think it's working, my confusion is caused because the output filenames don't match those in the tutorial example, I get fewer frames and an error window appears:

in the tutorial, the final files are named: final_ionized.pdb final_trajectory.dcd initial_ionized.pdb initial_trajectory.pdb

whereas what I get is: walker1_ionized.pdb walker1_trajectory.dcd walker2_ionized.pdb walker2_trajectory.dcd

they have just 6 frames in the dcd files rather than 41 (in the tutorial files)

the error window shows this (I used 'test' as the 'output prefix' leading to 'test.tcl'):


while executing

"exit" (file "/home/md/comd_tutorial_files/test.tcl" line 515) invoked from within "source $tcl_file_name" (procedure "::comd::Prepare_system" line 966) invoked from within "::comd::Prepare_system" invoked from within ".comdgui.main_frame.prepare.button invoke" ("uplevel" body line 1) invoked from within "uplevel #0 [list $w invoke]" (procedure "tk::ButtonUp" line 22) invoked from within "tk::ButtonUp .comdgui.main_frame.prepare.button" (command bound to event)


The parameters were:

Initial PDB 1ake - chain A Final PDB 4ake - chain A Box padding X 15A, Y 15A, Z 15A add counter ions - True Minimization 298K, 500 steps Step cutoff 4.0A Acceptance ratio 0.1 Deviation 0.1 Max ANM steps 100000 Spring constant 20000 TMD length 10ps coMD cycles 100

is it simply that using an acceptance ratio of 0.1 rather than 0.9 (which would allow for a more convoluted pathway) leads to very rapid convergence of the two structures ? (I was expecting everything to take much longer !)

BTW, I also had to delete the 'AP5' residue from 1ake.pdb to get things to run - is that normal ?

jamesmkrieger commented 5 years ago

The number of frames depends on the parameters used. I definitely need to update that tutorial and describe things more clearly. I get 6 cycles too.

As for the filenames, I changed those at some point because I was thinking there could be multiple walkers starting from the same structure. I don't think I ended up implementing it like that though.

I also got that error message but everything seemed to be fine so I'm not sure what's up with that.

A step cutoff of 4 A is way too big. I'd suggest one around 1 A for a protein of that size. Reducing that would definitely give you more steps. The other parameters seem good.

I think an acceptance ratio of 0.1 is ok. Yes, I think 0.9 would be make longer paths but isn't really necessary.

Thanks for the detailed feedback as well as the bug report.

Best wishes James

On Nov 28, 2019, at 5:41 PM, abvl notifications@github.com wrote:

 Dear James,

thanks for responding so quickly...

on reflection, I think it's working, my confusion is caused because the output filenames don't match those in the tutorial example, I get fewer frames and an error window appears:

in the tutorial, the final files are named: final_ionized.pdb final_trajectory.dcd initial_ionized.pdb initial_trajectory.pdb

whereas what I get is: walker1_ionized.pdb walker1_trajectory.dcd walker2_ionized.pdb walker2_trajectory.dcd

they have just 6 frames in the dcd files rather than 41 (in the tutorial files)

the error window shows this (I used 'test' as the 'output prefix' leading to 'test.tcl')

while executing "exit" (file "/home/md/comd_tutorial_files/test.tcl" line 515) invoked from within "source $tcl_file_name" (procedure "::comd::Prepare_system" line 966) invoked from within "::comd::Prepare_system" invoked from within ".comdgui.main_frame.prepare.button invoke" ("uplevel" body line 1) invoked from within "uplevel #0 [list $w invoke]" (procedure "tk::ButtonUp" line 22) invoked from within "tk::ButtonUp .comdgui.main_frame.prepare.button" (command bound to event)

The parameters were:

Initial PDB 1ake - chain A Final PDB 4ake - chain A Box padding X 15A, Y 15A, Z 15A add counter ions - True Minimization 298K, 500 steps Step cutoff 4.0A Acceptance ratio 0.1 Deviation 0.1 Max ANM steps 100000 Spring constant 20000 TMD length 10ps coMD cycles 100

is it simply that using an acceptance ratio of 0.1 rather than 0.9 (which would allow for a more convoluted pathway) leads to very rapid convergence of the two structures ? (I was expecting everything to take much longer !)

BTW, I also had to delete the 'AP5' residue from 1ake.pdb to get things to run - is that normal ?

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub, or unsubscribe.

jamesmkrieger commented 5 years ago

You probably don't need to exclude the ligand AP5 residue. It should be ignored anyway in the ANM part and there should be a way of accommodating it in the all-atom part. Did you include the Charmm general forcefield CGenFF?

abvl commented 5 years ago

Thanks for all your help with this, I now have everything working.

jamesmkrieger commented 5 years ago

You're welcome