[gmx-users] RMSD
Hi I have 2 questions about rmsd calculation: 1) why rmsd calculation is done on heavy atoms? 2) what means of mass weighted superposition (rmsd is calculated after mass weighted superposition) Any help will highly appreciated! -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] 6LYZ.pdb + Gromacs version 4.0.5 = Simulation Broken
Hi Lin, First of all, I would suggest sticking to a single processor until you have a protocol that works. Previously you had an issue with the addition of ions to your .top file. In your protocol, it's not mentioned. Have you made sure that issue is cleared? Cheers, Tsjerk On Sun, Jan 10, 2010 at 5:25 AM, Justin A. Lemkul jalem...@vt.edu wrote: On 1/9/10 10:42 PM, Chih-Ying Lin wrote: Hi I did the EM and the potential energy went to the very negative number. But the simulaiton broke in the PR step. Thank you Lin Saying the system broke is useless. There will certainly be some information in the .log file (LINCS warnings, etc). If that is the case, then you need to follow the standard advice that is always given in such cases: http://www.gromacs.org/Documentation/Errors#LINCS.2fSETTLE.2fSHAKE_warnings http://www.gromacs.org/Documentation/Terminology/Blowing_Up If you want more detailed advice, provide real output - information from how well the EM converged, messages in the .log file (unless it's just LINCS stuff, see the above links). Also realize that EM does not always remove all potentially problematic contacts. Your system may require a bit more finesse. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing list gmx-us...@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- Tsjerk A. Wassenaar, Ph.D. Computational Chemist Medicinal Chemist Neuropharmacologist -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] InflateGRO and trimer simulation
Hello Gromacs Users, I would like to run a simulation of a trimer in a DPPC membrane. I really like the elegant solution that inflategro script offers, however I'm afraid I won't be able to use it, because I need to have lipids in a small space between the monomers right in the center of the box and these will certainly get deleted during the whole procedure. I also need to keep the waters that come with the structure of the protein. Do you have any idea how I can elegantly embed my trimer in a membrane, possibly using somehow modified inflategro? Christopher -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] time-average structure
Hi what command is useful for obtaining time-average structure? Any help will highly appreciated! -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] interaction energy
Hi I want to obtain interaction energy between protein and dna in simulation pr-dna complex. what command is suitable for that? Any help will highly appreciated! -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] time-average structure
On 1/10/10 6:56 AM, leila karami wrote: Hi what command is useful for obtaining time-average structure? Something like g_cluster might do the trick, but be advised: http://www.gromacs.org/Documentation/Terminology/Average_Structure -Justin Any help will highly appreciated! -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] InflateGRO and trimer simulation
On 1/10/10 5:17 AM, KM wrote: Hello Gromacs Users, I would like to run a simulation of a trimer in a DPPC membrane. I really like the elegant solution that inflategro script offers, however I'm afraid I won't be able to use it, because I need to have lipids in a small space between the monomers right in the center of the box and these will certainly get deleted during the whole procedure. I also need to keep the waters that come with the structure of the protein. Do you have any idea how I can elegantly embed my trimer in a membrane, possibly using somehow modified inflategro? Christopher I've built multimeric systems with InflateGRO, but with a lot less constraints than you have :) You can, for instance, run an inflation step with a very minimal scaling factor (like 1.01 or 1.05) and a very small cutoff (7 or so, but you'd have to play around with this) so you essentially delete lipids in place without moving them too much. I don't know if the necessary lipids between the monomers will be affected. As for preserving water molecules, you could probably just extract their coordinates from a suitably-oriented starting structure and paste them into the InflateGRO output. None of this requires modifying InflateGRO, but will likely require a lot of trial and error, if it even works, given the very specific nature of what you need. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] RMSD
leila karami wrote: Hi I have 2 questions about rmsd calculation: 1) why rmsd calculation is done on heavy atoms? 2) what means of mass weighted superposition (rmsd is calculated after mass weighted superposition) Consider the difference between centre of mass and center of geometry... (Wikipedia is probably your friend here!) Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] InflateGRO and trimer simulation
Try the doughnut mode in the newest version of the program -- it's meant for exactly this situation. http://www.csb.bit.uni-bonn.de/inflategro.html -- original message -- Hello Gromacs Users, I would like to run a simulation of a trimer in a DPPC membrane. I really like the elegant solution that inflategro script offers, however I'm afraid I won't be able to use it, because I need to have lipids in a small space between the monomers right in the center of the box and these will certainly get deleted during the whole procedure. I also need to keep the waters that come with the structure of the protein. Do you have any idea how I can elegantly embed my trimer in a membrane, possibly using somehow modified inflategro? Christopher -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] InflateGRO and trimer simulation
Hi, the new inflategro with the 'doughnut mode' (http://www.csb.bit.uni-bonn.de/inflategro.html) might do what you want. Ciao, Patrick Justin A. Lemkul a écrit : On 1/10/10 5:17 AM, KM wrote: Hello Gromacs Users, I would like to run a simulation of a trimer in a DPPC membrane. I really like the elegant solution that inflategro script offers, however I'm afraid I won't be able to use it, because I need to have lipids in a small space between the monomers right in the center of the box and these will certainly get deleted during the whole procedure. I also need to keep the waters that come with the structure of the protein. Do you have any idea how I can elegantly embed my trimer in a membrane, possibly using somehow modified inflategro? Christopher I've built multimeric systems with InflateGRO, but with a lot less constraints than you have :) You can, for instance, run an inflation step with a very minimal scaling factor (like 1.01 or 1.05) and a very small cutoff (7 or so, but you'd have to play around with this) so you essentially delete lipids in place without moving them too much. I don't know if the necessary lipids between the monomers will be affected. As for preserving water molecules, you could probably just extract their coordinates from a suitably-oriented starting structure and paste them into the InflateGRO output. None of this requires modifying InflateGRO, but will likely require a lot of trial and error, if it even works, given the very specific nature of what you need. -Justin -- ___ new E-mail address: patrick.fu...@univ-paris-diderot.fr Patrick FUCHS Dynamique des Structures et Interactions des Macromolécules Biologiques INTS, INSERM UMR-S665, Université Paris Diderot, 6 rue Alexandre Cabanel, 75015 Paris Tel : +33 (0)1-44-49-30-57 - Fax : +33 (0)1-47-34-74-31 Web Site: http://www.dsimb.inserm.fr/~fuchs -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Exceeding of Maximum allowed number of DD cells
Dear GMX-Users, I'm testing my 256 full hydrated lipid on blue gene. The purpose is to find out the right number for -npme, as mdrun can not estimate itself successfully. I met the problem that how to match the maxinum allowed number for DD cells with large number of CPU cores. My simulation box size is about 8x8x9nm^3, with normal LINCS parameter and dds=0.8. The log file said that the maximum allowed number of DD cells is 8x8x9. As far as I understand, DD assigns one core to one cell, so the maximum core I can use in this case for PP part is 8x8x9=576 cores. I then ran with 512 cores with -npme=128. My system runs without problem. What if I want to use more cores? Then I try to increase the -dds from 0.8 to 0.9, this leads to an increasing of the maximum allowed number of DD cells to 8x10x10. This time is 1024 cores in total and I set -npme=224, then PP part will have 800 cores which are within 10x10x10. The system ran initially but corrupted very soon with warning that DD cell 2 1 4 could on obtain 56 of the 57 atoms that are connected via constraints from the neighboring cells Therefore the dilemma is if I increase the -dds, I can meet the requirement for the maximum allowed number of DD cells, but fail the maximum length of constraints in LINCS. Does it mean that for a relative small system, it is not possible to using up to thousand of cores by domain decomposition? I know that if it makes more sense to use thousand of cores for huge system, but if my purpose is simply to speed up the simulation, what should I do? Thank you. Chao -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] free energy
Hello Justin, Should I add dum_opls no. in atom type file in which Lennard-jones parameters will set to zero. OR. I don't have to modify anything prog.will take care of lennard-jones interactions. I have put charge zero on all atom in TYPE B. Nilesh On Sun, January 3, 2010 12:57 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: THanks Justin. I am using Groamcs 4.0.7 version. I will do more simulation with more lambada value betwen 0 to 0.1. I really think this is less significant than point #1 I mentioned before. Be sure you are actually setting the relevant .mdp parameters (couple-lambda0, couple-lambda1, etc) as they will heavily influence the results. Leaving them as default (as you are doing, per the .mdp file you provided) may give you something other than what you want. And do consult the links I provided; they give excellent background on the nature of these calculations. -Justin Thanks. Nilesh On Sun, January 3, 2010 12:39 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: Hey All, I am trying to calculate the solvation free energy of glucose in ionic liquids. I am getting really large value for lamda=0 in vacuum and solvent as well. I am getting final solvation energy value ~ 310 kJ/mol. Can you tell me where I am wrong. There are several possibilities: 1. Are you decoupling LJ and Coulombic interactions simultaneously? They should be done in separate steps, or else you get very large errors. 2. You probably need a lot more lambda points to properly define the curve. For strongly hydrogen-bonding solutes, you should probably include closely-spaced points near lambda = 0 (the fully-coupled state), like 0.005, 0.01, 0.015, etc. The curve from 0 to 0.1 will otherwise be ill-defined and prone to large errors. You also do not say which version of the code you are using. The free energy code has been tweaked in version 4.0 such that topology handling is different, and the procedure is actually quite a bit more streamlined, given the .mdp options. Might I recommend: http://www.alchemistry.org/wiki/index.php/Best_Practices http://md.chem.rug.nl/education/Free-Energy_Course/index.html -Justin Here I have pasted the the lamda value in vacuum and solvent Vacuum Solvent 0.0 527.16 0.0 4284.769 0.1 74.119 0.1 161.606 0.2 37.96650.2 -16.3176 0.3 28.748 0.3 -52.488 0.4 7.555 0.4 -4.3341 0.5 1.495 0.5 6.996 0.6 -4.796 0.6 155.485 0.7 2.361 0.7 174.20 0.8 3.69 0.8 193.016 0.9 5.483 0.9 201.820 0.92 6.391 0.92 211.79 1.0 8.887 1.0 227.79 For each lamda value I equibrate the system for 1ns and run the dynamics for 4 ns. here is dynamics file I used. pasted the topology file at the end. title = cpeptide MD cpp = /usr/bin/cpp constraints = none integrator = md dt = 0.001 ; ps ! nsteps = 400 ; total 5 ps. nstcomm = 1 nstxout = 50 nstvout = 0 nstfout = 0 nstlist = 10 ns_type = grid rlist = 1.0 rcoulomb= 1.0 rvdw= 1.4 coulombtype = PME vdwtype = cut-off pbc = xyz fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes ; Berendsen temperature coupling is on Tcoupl = v-rescale tau_t = 0.1 tc-grps =system ref_t = 350 ; Pressure coupling is on Pcoupl = parrinello-rahman pcoupltype = isotropic tau_p = 0.5 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocites is on at 300 K. gen_vel = yes gen_temp= 350.0 gen_seed = 173529 free_energy = yes init_lambda = 0.00 sc-alpha = 0.6 sc_power= 1 Thanks Nilesh Topology file (not pasted whole because of size) ; ; File 'glu-ac.top' was generated ; By user: ndhumal (36026) ; On host: c20 ; At date: Tue Dec 15 15:25:07 2009 ; ; This is your topology file ; Gromacs Runs One Microsecond At Cannonball Speeds ; ; Include forcefield parameters #include ffoplsaa.itp [ moleculetype ] ; Namenrexcl Protein 3 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_180 1GUL O 1 -0.415.9994 opls_180 0.015.9994 2 opls_159 1GUL C2 2 0.365 12.011 opls_159 0.0 12.011 3 opls_159 1GUL C3 2 0.205 12.011 opls_159 0.0 12.011 4 opls_159 1GUL C4 2 0.205 12.011 opls_159 0.0 12.011 5 opls_159 1GUL C5 2 0.205 12.011 opls_159 0.0 12.011 6 opls_137 1GUL C6 2 0.17 12.011 opls_137 0.0 12.011 7 opls_159 1GUL C7 2 0.145 12.011 opls_159 0.0 12.011
Re: [gmx-users] free energy
On 1/10/10 1:03 PM, Nilesh Dhumal wrote: Hello Justin, Should I add dum_opls no. in atom type file in which Lennard-jones parameters will set to zero. OR. I don't have to modify anything prog.will take care of lennard-jones interactions. I have put charge zero on all atom in TYPE B. You shouldn't have to do make any modifications to the topology whatsoever, according to the documentation. The couple-lambda0 and couple-lambda1 parameters should do the decoupling for you. -Justin Nilesh On Sun, January 3, 2010 12:57 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: THanks Justin. I am using Groamcs 4.0.7 version. I will do more simulation with more lambada value betwen 0 to 0.1. I really think this is less significant than point #1 I mentioned before. Be sure you are actually setting the relevant .mdp parameters (couple-lambda0, couple-lambda1, etc) as they will heavily influence the results. Leaving them as default (as you are doing, per the .mdp file you provided) may give you something other than what you want. And do consult the links I provided; they give excellent background on the nature of these calculations. -Justin Thanks. Nilesh On Sun, January 3, 2010 12:39 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: Hey All, I am trying to calculate the solvation free energy of glucose in ionic liquids. I am getting really large value for lamda=0 in vacuum and solvent as well. I am getting final solvation energy value ~ 310 kJ/mol. Can you tell me where I am wrong. There are several possibilities: 1. Are you decoupling LJ and Coulombic interactions simultaneously? They should be done in separate steps, or else you get very large errors. 2. You probably need a lot more lambda points to properly define the curve. For strongly hydrogen-bonding solutes, you should probably include closely-spaced points near lambda = 0 (the fully-coupled state), like 0.005, 0.01, 0.015, etc. The curve from 0 to 0.1 will otherwise be ill-defined and prone to large errors. You also do not say which version of the code you are using. The free energy code has been tweaked in version 4.0 such that topology handling is different, and the procedure is actually quite a bit more streamlined, given the .mdp options. Might I recommend: http://www.alchemistry.org/wiki/index.php/Best_Practices http://md.chem.rug.nl/education/Free-Energy_Course/index.html -Justin Here I have pasted the the lamda value in vacuum and solvent Vacuum Solvent 0.0 527.16 0.0 4284.769 0.1 74.119 0.1 161.606 0.2 37.9665 0.2 -16.3176 0.3 28.748 0.3 -52.488 0.4 7.555 0.4 -4.3341 0.5 1.495 0.5 6.996 0.6 -4.796 0.6 155.485 0.7 2.361 0.7 174.20 0.8 3.690.8 193.016 0.9 5.483 0.9 201.820 0.92 6.391 0.92 211.79 1.0 8.887 1.0 227.79 For each lamda value I equibrate the system for 1ns and run the dynamics for 4 ns. here is dynamics file I used. pasted the topology file at the end. title = cpeptide MD cpp = /usr/bin/cpp constraints = none integrator = md dt = 0.001 ; ps ! nsteps = 400 ; total 5 ps. nstcomm = 1 nstxout = 50 nstvout = 0 nstfout = 0 nstlist = 10 ns_type = grid rlist = 1.0 rcoulomb= 1.0 rvdw= 1.4 coulombtype = PME vdwtype = cut-off pbc = xyz fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes ; Berendsen temperature coupling is on Tcoupl = v-rescale tau_t = 0.1 tc-grps =system ref_t = 350 ; Pressure coupling is on Pcoupl = parrinello-rahman pcoupltype = isotropic tau_p = 0.5 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocites is on at 300 K. gen_vel = yes gen_temp= 350.0 gen_seed = 173529 free_energy = yes init_lambda = 0.00 sc-alpha = 0.6 sc_power= 1 Thanks Nilesh Topology file (not pasted whole because of size) ; ; File 'glu-ac.top' was generated ; By user: ndhumal (36026) ; On host: c20 ; At date: Tue Dec 15 15:25:07 2009 ; ; This is your topology file ; Gromacs Runs One Microsecond At Cannonball Speeds ; ; Include forcefield parameters #include ffoplsaa.itp [ moleculetype ] ; Namenrexcl Protein 3 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_180 1GUL O 1 -0.415.9994 opls_180 0.015.9994 2 opls_159 1GUL C2 2 0.365 12.011 opls_159 0.0 12.011 3 opls_159 1GUL C3 2 0.205 12.011 opls_159 0.0 12.011 4 opls_159 1GUL C4 2 0.205 12.011 opls_159 0.0 12.011 5 opls_159 1
Re: [gmx-users] InflateGRO and trimer simulation
Thank you! This is what I need, however I keep getting a lot of errors about uninitialized values. When the script finishes, the membrane is rescaled, but the protein is untouched in the corner of the new box :( I'll try to investigate a little bit despite I don't know perl, however I believe that knowledge of python, C and tcl will help me understand what's going on. Chris 2010/1/10 patrick fuchs patrick.fu...@univ-paris-diderot.fr: Hi, the new inflategro with the 'doughnut mode' (http://www.csb.bit.uni-bonn.de/inflategro.html) might do what you want. Ciao, Patrick Justin A. Lemkul a écrit : On 1/10/10 5:17 AM, KM wrote: Hello Gromacs Users, I would like to run a simulation of a trimer in a DPPC membrane. I really like the elegant solution that inflategro script offers, however I'm afraid I won't be able to use it, because I need to have lipids in a small space between the monomers right in the center of the box and these will certainly get deleted during the whole procedure. I also need to keep the waters that come with the structure of the protein. Do you have any idea how I can elegantly embed my trimer in a membrane, possibly using somehow modified inflategro? Christopher I've built multimeric systems with InflateGRO, but with a lot less constraints than you have :) You can, for instance, run an inflation step with a very minimal scaling factor (like 1.01 or 1.05) and a very small cutoff (7 or so, but you'd have to play around with this) so you essentially delete lipids in place without moving them too much. I don't know if the necessary lipids between the monomers will be affected. As for preserving water molecules, you could probably just extract their coordinates from a suitably-oriented starting structure and paste them into the InflateGRO output. None of this requires modifying InflateGRO, but will likely require a lot of trial and error, if it even works, given the very specific nature of what you need. -Justin -- ___ new E-mail address: patrick.fu...@univ-paris-diderot.fr Patrick FUCHS Dynamique des Structures et Interactions des Macromolécules Biologiques INTS, INSERM UMR-S665, Université Paris Diderot, 6 rue Alexandre Cabanel, 75015 Paris Tel : +33 (0)1-44-49-30-57 - Fax : +33 (0)1-47-34-74-31 Web Site: http://www.dsimb.inserm.fr/~fuchs -- gmx-users mailing list gmx-us...@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- W każdej wsi jest pochodnia oświaty - nauczyciel - oraz gaśnica - ksiądz. Victor Hugo -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] InflateGRO and trimer simulation
OK, now it works. :) Chris W dniu 10 stycznia 2010 20:06 użytkownik KM mitomas...@gmail.com napisał: Thank you! This is what I need, however I keep getting a lot of errors about uninitialized values. When the script finishes, the membrane is rescaled, but the protein is untouched in the corner of the new box :( I'll try to investigate a little bit despite I don't know perl, however I believe that knowledge of python, C and tcl will help me understand what's going on. Chris 2010/1/10 patrick fuchs patrick.fu...@univ-paris-diderot.fr: Hi, the new inflategro with the 'doughnut mode' (http://www.csb.bit.uni-bonn.de/inflategro.html) might do what you want. Ciao, Patrick Justin A. Lemkul a écrit : On 1/10/10 5:17 AM, KM wrote: Hello Gromacs Users, I would like to run a simulation of a trimer in a DPPC membrane. I really like the elegant solution that inflategro script offers, however I'm afraid I won't be able to use it, because I need to have lipids in a small space between the monomers right in the center of the box and these will certainly get deleted during the whole procedure. I also need to keep the waters that come with the structure of the protein. Do you have any idea how I can elegantly embed my trimer in a membrane, possibly using somehow modified inflategro? Christopher I've built multimeric systems with InflateGRO, but with a lot less constraints than you have :) You can, for instance, run an inflation step with a very minimal scaling factor (like 1.01 or 1.05) and a very small cutoff (7 or so, but you'd have to play around with this) so you essentially delete lipids in place without moving them too much. I don't know if the necessary lipids between the monomers will be affected. As for preserving water molecules, you could probably just extract their coordinates from a suitably-oriented starting structure and paste them into the InflateGRO output. None of this requires modifying InflateGRO, but will likely require a lot of trial and error, if it even works, given the very specific nature of what you need. -Justin -- ___ new E-mail address: patrick.fu...@univ-paris-diderot.fr Patrick FUCHS Dynamique des Structures et Interactions des Macromolécules Biologiques INTS, INSERM UMR-S665, Université Paris Diderot, 6 rue Alexandre Cabanel, 75015 Paris Tel : +33 (0)1-44-49-30-57 - Fax : +33 (0)1-47-34-74-31 Web Site: http://www.dsimb.inserm.fr/~fuchs -- gmx-users mailing list gmx-us...@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- W każdej wsi jest pochodnia oświaty - nauczyciel - oraz gaśnica - ksiądz. Victor Hugo -- W każdej wsi jest pochodnia oświaty - nauczyciel - oraz gaśnica - ksiądz. Victor Hugo -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Unstable Minimizations
I am trying to get this workflow opperational. However, my systems are getting unstable. I have preped two mdp files: 1) one for restrained 2) unrestrained. LINCS errors appear for restrained and unrestrained has infinite energy appearing. http://boinc.drugdiscoveryathome.com/*em_restrained_rcs_mdrun.txthttp://boinc.drugdiscoveryathome.com/em_restrained_rcs_mdrun.txt * ** This is where I get the LINCS Warnings Step -1, time -0.001 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.461520, max 14.428611 (between atoms 1668 and 1669) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length Steepest Descents: Tolerance (Fmax) = 1.0e+04 Number of steps= 100 Warning: 1-4 interaction between 1658 and 1672 at distance 2.655 which is larger than the 1-4 table size 2.400 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size Step=0, Dmax= 1.0e-02 nm, Epot= 1.09364e+09 Fmax= 2.21154e+11, atom= 3292 Step 1, time 0.001 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.680509, max 22.293625 (between atoms 1668 and 1670) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length What is a reasonable increase in table-extension. Is this a mis-leading suggestion? Here is the log from the unrestrained minimization. http://boinc.drugdiscoveryathome.com/*em_rcs_mdrun.txthttp://boinc.drugdiscoveryathome.com/em_rcs_mdrun.txt * Here is a zip archive containing the working directory for this minimization. Its about 428 kb http://boinc.drugdiscoveryathome.com/rcs_ga_run_10_bt_Fzd2-MD7-MD8-7.zip_lig_24205_ChemDiv_5754-2873_ts_1263004110202172000.zip -- Jack http://drugdiscoveryathome.com http://hydrogenathome.org -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Unstable Minimizations
On 1/10/10 5:18 PM, Jack Shultz wrote: I am trying to get this workflow opperational. However, my systems are getting unstable. I have preped two mdp files: 1) one for restrained 2) unrestrained. LINCS errors appear for restrained and unrestrained has infinite energy appearing. http://boinc.drugdiscoveryathome.com/_em_restrained_rcs_mdrun.txt_ http://boinc.drugdiscoveryathome.com/em_restrained_rcs_mdrun.txt __ This log file shows several long bond warnings, which may be the root of your problem. See here: http://www.gromacs.org/Documentation/Errors#Long_bonds_and.2for_missing_atoms Since your minimization is failing immediately, there is something physically unreasonable about your structure, such that EM cannot resolve the problem. Note, too, that one of the long bond warnings pertained to atom 1668, which is the location of the first LINCS warning. Coincidence? Not likely. Re-examine the starting structure and figure out if anything is missing or poorly reconstructed (e.g., from initially missing atoms). This is where I get the LINCS Warnings Step -1, time -0.001 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.461520, max 14.428611 (between atoms 1668 and 1669) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length Steepest Descents: Tolerance (Fmax) = 1.0e+04 Number of steps= 100 Warning: 1-4 interaction between 1658 and 1672 at distance 2.655 which is larger than the 1-4 table size 2.400 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size Step=0, Dmax= 1.0e-02 nm, Epot= 1.09364e+09 Fmax= 2.21154e+11, atom= 3292 Step 1, time 0.001 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.680509, max 22.293625 (between atoms 1668 and 1670) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length What is a reasonable increase in table-extension. Is this a mis-leading suggestion? You should not adjust the table-extension. The other part of the error message is what you need to pay attention to (your system is exploding). -Justin Here is the log from the unrestrained minimization. http://boinc.drugdiscoveryathome.com/_em_rcs_mdrun.txt_ http://boinc.drugdiscoveryathome.com/em_rcs_mdrun.txt Here is a zip archive containing the working directory for this minimization. Its about 428 kb http://boinc.drugdiscoveryathome.com/rcs_ga_run_10_bt_Fzd2-MD7-MD8-7.zip_lig_24205_ChemDiv_5754-2873_ts_1263004110202172000.zip -- Jack http://drugdiscoveryathome.com http://hydrogenathome.org -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Retinal and Schiff base parameters for ffG53a6
Hello, I would like to simulate bacteriorhodopsin embedded into a membrane. I think ffG53a6 is a good choice - force field is quite new (at least newer than ffgmx) and there are topologies available for the most popular lipids, for example the set developed by Andreas Kukol. To run a simulation I will need parameters for retinal and Schiff base connecting Lys216 with retinal. I have those parameters for ffgmx, however they are incompatible with ffG53a6 (and that's no surprise). Does anyone know where I can get parameters for ffG53a6 force field (or compatible with it)? Chris -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Retinal and Schiff base parameters for ffG53a6
On 1/10/10 6:16 PM, KM wrote: Hello, I would like to simulate bacteriorhodopsin embedded into a membrane. I think ffG53a6 is a good choice - force field is quite new (at least newer than ffgmx) and there are topologies available for the most popular lipids, for example the set developed by Andreas Kukol. To run a simulation I will need parameters for retinal and Schiff base connecting Lys216 with retinal. I have those parameters for ffgmx, however they are incompatible with ffG53a6 (and that's no surprise). Does anyone know where I can get parameters for ffG53a6 force field (or compatible with it)? Chris If they're not already in the literature or the User Contribution site on the Gromacs website, you'll likely have to develop them yourself. See, for example: http://www.gromacs.org/Documentation/How-tos/Parameterization I know I have seen several simulations of backteriorhodopsin in the literature, so perhaps you can find parameters that have already been published. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] free energy
Thanks Justin, I am not setting any Lennard-Jones parameters to zero for part B. I put zero charge on all atoms for Part B. (Lennard Jones parameters for B are same as Part A). I am using default parameters for couple-lambda0 and couple-lambda1. For one of my simulation (glucose + ionic liquids) I am getting solvation free energy 35.5 Kcal/mol. Is it too high? Nilesh On Sun, January 10, 2010 1:08 pm, Justin A. Lemkul wrote: On 1/10/10 1:03 PM, Nilesh Dhumal wrote: Hello Justin, Should I add dum_opls no. in atom type file in which Lennard-jones parameters will set to zero. OR. I don't have to modify anything prog.will take care of lennard-jones interactions. I have put charge zero on all atom in TYPE B. You shouldn't have to do make any modifications to the topology whatsoever, according to the documentation. The couple-lambda0 and couple-lambda1 parameters should do the decoupling for you. -Justin Nilesh On Sun, January 3, 2010 12:57 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: THanks Justin. I am using Groamcs 4.0.7 version. I will do more simulation with more lambada value betwen 0 to 0.1. I really think this is less significant than point #1 I mentioned before. Be sure you are actually setting the relevant .mdp parameters (couple-lambda0, couple-lambda1, etc) as they will heavily influence the results. Leaving them as default (as you are doing, per the .mdp file you provided) may give you something other than what you want. And do consult the links I provided; they give excellent background on the nature of these calculations. -Justin Thanks. Nilesh On Sun, January 3, 2010 12:39 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: Hey All, I am trying to calculate the solvation free energy of glucose in ionic liquids. I am getting really large value for lamda=0 in vacuum and solvent as well. I am getting final solvation energy value ~ 310 kJ/mol. Can you tell me where I am wrong. There are several possibilities: 1. Are you decoupling LJ and Coulombic interactions simultaneously? They should be done in separate steps, or else you get very large errors. 2. You probably need a lot more lambda points to properly define the curve. For strongly hydrogen-bonding solutes, you should probably include closely-spaced points near lambda = 0 (the fully-coupled state), like 0.005, 0.01, 0.015, etc. The curve from 0 to 0.1 will otherwise be ill-defined and prone to large errors. You also do not say which version of the code you are using. The free energy code has been tweaked in version 4.0 such that topology handling is different, and the procedure is actually quite a bit more streamlined, given the .mdp options. Might I recommend: http://www.alchemistry.org/wiki/index.php/Best_Practices http://md.chem.rug.nl/education/Free-Energy_Course/index.html -Justin Here I have pasted the the lamda value in vacuum and solvent VacuumSolvent 0.0 527.16 0.0 4284.769 0.1 74.119 0.1 161.606 0.2 37.9665 0.2 -16.3176 0.3 28.748 0.3 -52.488 0.4 7.5550.4 -4.3341 0.5 1.4950.5 6.996 0.6 -4.796 0.6 155.485 0.7 2.3610.7 174.20 0.8 3.69 0.8 193.016 0.9 5.4830.9 201.820 0.92 6.391 0.92 211.79 1.0 8.8871.0 227.79 For each lamda value I equibrate the system for 1ns and run the dynamics for 4 ns. here is dynamics file I used. pasted the topology file at the end. title = cpeptide MD cpp = /usr/bin/cpp constraints = none integrator = md dt = 0.001 ; ps ! nsteps = 400 ; total 5 ps. nstcomm = 1 nstxout = 50 nstvout = 0 nstfout = 0 nstlist = 10 ns_type = grid rlist = 1.0 rcoulomb= 1.0 rvdw= 1.4 coulombtype = PME vdwtype = cut-off pbc = xyz fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes ; Berendsen temperature coupling is on Tcoupl = v-rescale tau_t = 0.1 tc-grps =system ref_t = 350 ; Pressure coupling is on Pcoupl = parrinello-rahman pcoupltype = isotropic tau_p = 0.5 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocites is on at 300 K. gen_vel = yes gen_temp = 350.0 gen_seed = 173529 free_energy = yes init_lambda = 0.00 sc-alpha = 0.6 sc_power= 1 Thanks Nilesh Topology file (not pasted whole because of size) ; ;File 'glu-ac.top' was generated ;By user: ndhumal (36026) ;On host: c20 ;At date: Tue Dec 15 15:25:07 2009 ; ;This is your topology file ;Gromacs Runs One Microsecond At Cannonball Speeds ; ; Include forcefield parameters
Re: [gmx-users] Exceeding of Maximum allowed number of DD cells
Chao Zhang wrote: Dear GMX-Users, I'm testing my 256 full hydrated lipid on blue gene. The purpose is to find out the right number for -npme, as mdrun can not estimate itself successfully. I met the problem that how to match the maxinum allowed number for DD cells with large number of CPU cores. My simulation box size is about 8x8x9nm^3, with normal LINCS parameter and dds=0.8. The log file said that the maximum allowed number of DD cells is 8x8x9. Why are you setting -dds? -rcon and -rdd are also variables to play with... but if you get too close to the bone you can run into (e.g.) LINCS problems. A lipid-water system has inhomogeneous interaction density, so the DD load-balance needs to scale the starting guess for cells, and -dds will have a significant effect at high parallelization. See mdrun -h, and the manual and GROMACS 4 paper. As far as I understand, DD assigns one core to one cell, so the maximum core I can use in this case for PP part is 8x8x9=576 cores. I then ran with 512 cores with -npme=128. My system runs without problem. What if I want to use more cores? Then I try to increase the -dds from 0.8 to 0.9, this leads to an increasing of the maximum allowed number of DD cells to 8x10x10. This time is 1024 cores in total and I set -npme=224, then PP part will have 800 cores which are within 10x10x10. The system ran initially but corrupted very soon with warning that DD cell 2 1 4 could on obtain 56 of the 57 atoms that are connected via constraints from the neighboring cells Therefore the dilemma is if I increase the -dds, I can meet the requirement for the maximum allowed number of DD cells, but fail the maximum length of constraints in LINCS. Does it mean that for a relative small system, it is not possible to using up to thousand of cores by domain decomposition? Yes. Each cell takes responsibility for a subset of atoms, and then communicates them to neighbouring cell who need to know. As the cell gets smaller, the communication cost would get larger. GROMACS sets a number of semi-artificial constraints on the cell size with the above options. There is a lower limit on DD cell size in practice for a given system with the GROMACS 4 implementation, but you have to experiment to find it. Whether you derive any speed advantage from moving towards that limit will depend on the relative performance of your processors and network. IBM's Blue Matter MD code is supposed to work down at around 1 atom per core, but GROMACS isn't built to do that. I know that if it makes more sense to use thousand of cores for huge system, but if my purpose is simply to speed up the simulation, what should I do? If your objective is increasing effective sampling, REMD of 16 replicas of 64 cores (or similar) makes much sense. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] free energy
On 1/10/10 7:02 PM, Nilesh Dhumal wrote: Thanks Justin, I am not setting any Lennard-Jones parameters to zero for part B. I put zero charge on all atoms for Part B. (Lennard Jones parameters for B are same as Part A). Under Gromacs 4.0.x, you should be able to leave the topology alone entirely. No B state should need to be defined explicitly. I am using default parameters for couple-lambda0 and couple-lambda1. The default behavior is vdw-q, is it not? Then what you're doing still doesn't make sense to me entirely, because no decoupling would even be specified if that is the case. Can you post your actual .mdp file (mdout.mdp might be better if you're leaving it up to grompp to decide the defaults, etc)? I'll admit I haven't done much with the 4.0 free energy code, but it seems to me that with an unaltered topology (no B state explicitly defined) you should be able to set the following in the .mdp file: couple-lambda0 = vdw-q couple-lambda1 = vdw ...to decouple only Coulombic interactions, and likewise: couple-lambda0 = vdw couple-lambda1 = none ...to decouple van der Waals. Perhaps someone can correct me if I've gotten this wrong. For one of my simulation (glucose + ionic liquids) I am getting solvation free energy 35.5 Kcal/mol. Is it too high? Not a clue. Isn't that what you're trying to figure out? Is there any sort of literature value for such a parameter? In any case, it's more plausible than your original value. -Justin Nilesh On Sun, January 10, 2010 1:08 pm, Justin A. Lemkul wrote: On 1/10/10 1:03 PM, Nilesh Dhumal wrote: Hello Justin, Should I add dum_opls no. in atom type file in which Lennard-jones parameters will set to zero. OR. I don't have to modify anything prog.will take care of lennard-jones interactions. I have put charge zero on all atom in TYPE B. You shouldn't have to do make any modifications to the topology whatsoever, according to the documentation. The couple-lambda0 and couple-lambda1 parameters should do the decoupling for you. -Justin Nilesh On Sun, January 3, 2010 12:57 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: THanks Justin. I am using Groamcs 4.0.7 version. I will do more simulation with more lambada value betwen 0 to 0.1. I really think this is less significant than point #1 I mentioned before. Be sure you are actually setting the relevant .mdp parameters (couple-lambda0, couple-lambda1, etc) as they will heavily influence the results. Leaving them as default (as you are doing, per the .mdp file you provided) may give you something other than what you want. And do consult the links I provided; they give excellent background on the nature of these calculations. -Justin Thanks. Nilesh On Sun, January 3, 2010 12:39 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: Hey All, I am trying to calculate the solvation free energy of glucose in ionic liquids. I am getting really large value for lamda=0 in vacuum and solvent as well. I am getting final solvation energy value ~ 310 kJ/mol. Can you tell me where I am wrong. There are several possibilities: 1. Are you decoupling LJ and Coulombic interactions simultaneously? They should be done in separate steps, or else you get very large errors. 2. You probably need a lot more lambda points to properly define the curve. For strongly hydrogen-bonding solutes, you should probably include closely-spaced points near lambda = 0 (the fully-coupled state), like 0.005, 0.01, 0.015, etc. The curve from 0 to 0.1 will otherwise be ill-defined and prone to large errors. You also do not say which version of the code you are using. The free energy code has been tweaked in version 4.0 such that topology handling is different, and the procedure is actually quite a bit more streamlined, given the .mdp options. Might I recommend: http://www.alchemistry.org/wiki/index.php/Best_Practices http://md.chem.rug.nl/education/Free-Energy_Course/index.html -Justin Here I have pasted the the lamda value in vacuum and solvent Vacuum Solvent 0.0 527.16 0.0 4284.769 0.1 74.119 0.1 161.606 0.2 37.9665 0.2 -16.3176 0.3 28.748 0.3 -52.488 0.4 7.555 0.4 -4.3341 0.5 1.495 0.5 6.996 0.6 -4.796 0.6 155.485 0.7 2.361 0.7 174.20 0.8 3.690.8 193.016 0.9 5.483 0.9 201.820 0.92 6.391 0.92 211.79 1.0 8.887 1.0 227.79 For each lamda value I equibrate the system for 1ns and run the dynamics for 4 ns. here is dynamics file I used. pasted the topology file at the end. title = cpeptide MD cpp = /usr/bin/cpp constraints = none integrator = md dt = 0.001 ; ps ! nsteps = 400 ; total 5 ps. nstcomm = 1 nstxout = 50 nstvout = 0 nstfout = 0 nstlist = 10 ns_type = grid rlist = 1.0 rcoulomb= 1.0 rvdw= 1.4
Re: [gmx-users] free energy
Justin, Here I have pasted the data about free energy from mdout.mdp file ; Free energy control stuff free_energy = yes init_lambda = 0.20 delta-lambda = 0 sc-alpha = 0.6 sc_power = 1 sc-sigma = 0.3 couple-moltype = couple-lambda0 = vdw-q couple-lambda1 = vdw-q couple-intramol = no Nilesh On Sun, January 10, 2010 7:40 pm, Justin A. Lemkul wrote: On 1/10/10 7:02 PM, Nilesh Dhumal wrote: Thanks Justin, I am not setting any Lennard-Jones parameters to zero for part B. I put zero charge on all atoms for Part B. (Lennard Jones parameters for B are same as Part A). Under Gromacs 4.0.x, you should be able to leave the topology alone entirely. No B state should need to be defined explicitly. I am using default parameters for couple-lambda0 and couple-lambda1. The default behavior is vdw-q, is it not? Then what you're doing still doesn't make sense to me entirely, because no decoupling would even be specified if that is the case. Can you post your actual .mdp file (mdout.mdp might be better if you're leaving it up to grompp to decide the defaults, etc)? I'll admit I haven't done much with the 4.0 free energy code, but it seems to me that with an unaltered topology (no B state explicitly defined) you should be able to set the following in the .mdp file: couple-lambda0 = vdw-q couple-lambda1 = vdw ...to decouple only Coulombic interactions, and likewise: couple-lambda0 = vdw couple-lambda1 = none ...to decouple van der Waals. Perhaps someone can correct me if I've gotten this wrong. For one of my simulation (glucose + ionic liquids) I am getting solvation free energy 35.5 Kcal/mol. Is it too high? Not a clue. Isn't that what you're trying to figure out? Is there any sort of literature value for such a parameter? In any case, it's more plausible than your original value. -Justin Nilesh On Sun, January 10, 2010 1:08 pm, Justin A. Lemkul wrote: On 1/10/10 1:03 PM, Nilesh Dhumal wrote: Hello Justin, Should I add dum_opls no. in atom type file in which Lennard-jones parameters will set to zero. OR. I don't have to modify anything prog.will take care of lennard-jones interactions. I have put charge zero on all atom in TYPE B. You shouldn't have to do make any modifications to the topology whatsoever, according to the documentation. The couple-lambda0 and couple-lambda1 parameters should do the decoupling for you. -Justin Nilesh On Sun, January 3, 2010 12:57 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: THanks Justin. I am using Groamcs 4.0.7 version. I will do more simulation with more lambada value betwen 0 to 0.1. I really think this is less significant than point #1 I mentioned before. Be sure you are actually setting the relevant .mdp parameters (couple-lambda0, couple-lambda1, etc) as they will heavily influence the results. Leaving them as default (as you are doing, per the .mdp file you provided) may give you something other than what you want. And do consult the links I provided; they give excellent background on the nature of these calculations. -Justin Thanks. Nilesh On Sun, January 3, 2010 12:39 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: Hey All, I am trying to calculate the solvation free energy of glucose in ionic liquids. I am getting really large value for lamda=0 in vacuum and solvent as well. I am getting final solvation energy value ~ 310 kJ/mol. Can you tell me where I am wrong. There are several possibilities: 1. Are you decoupling LJ and Coulombic interactions simultaneously? They should be done in separate steps, or else you get very large errors. 2. You probably need a lot more lambda points to properly define the curve. For strongly hydrogen-bonding solutes, you should probably include closely-spaced points near lambda = 0 (the fully-coupled state), like 0.005, 0.01, 0.015, etc. The curve from 0 to 0.1 will otherwise be ill-defined and prone to large errors. You also do not say which version of the code you are using. The free energy code has been tweaked in version 4.0 such that topology handling is different, and the procedure is actually quite a bit more streamlined, given the .mdp options. Might I recommend: http://www.alchemistry.org/wiki/index.php/Best_Practices http://md.chem.rug.nl/education/Free-Energy_Course/index.html -Justin Here I have pasted the the lamda value in vacuum and solvent Vacuum Solvent 0.0 527.16 0.0 4284.769 0.1 74.119 0.1 161.606 0.2 37.96650.2 -16.3176 0.3 28.748 0.3 -52.488 0.4 7.555 0.4 -4.3341 0.5 1.495 0.5 6.996 0.6 -4.796 0.6 155.485 0.7 2.361 0.7 174.20 0.8 3.69 0.8 193.016 0.9 5.483 0.9 201.820 0.92 6.391 0.92 211.79 1.0 8.887 1.0
Re: [gmx-users] free energy
On 1/10/10 7:53 PM, Nilesh Dhumal wrote: Justin, Here I have pasted the data about free energy from mdout.mdp file ; Free energy control stuff free_energy = yes init_lambda = 0.20 delta-lambda = 0 sc-alpha = 0.6 sc_power = 1 sc-sigma = 0.3 couple-moltype = couple-lambda0 = vdw-q couple-lambda1 = vdw-q couple-intramol = no According the manual, you're not doing anything. Setting both couple-lambda0 and couple-lambda1 to vdw-q means you're not decoupling either van der Waals or Coulombic interactions. It may be working if the B state you've specifically defined contains zero charges, but this is a very roundabout (and potentially unreliable) way of doing things. Do you know what your answer corresponds to if you have a B state in the topology (which is the obsolete method as of version 4.0) but the .mdp file specifies no actual decoupling? Sounds unreliable to me. If anything, you may have managed to decouple the Coulombic terms only, but I still wouldn't trust it. Try the .mdp settings I posted before, with an unaltered topology. That would be the easiest way, the by the book method, so to speak. -Justin Nilesh On Sun, January 10, 2010 7:40 pm, Justin A. Lemkul wrote: On 1/10/10 7:02 PM, Nilesh Dhumal wrote: Thanks Justin, I am not setting any Lennard-Jones parameters to zero for part B. I put zero charge on all atoms for Part B. (Lennard Jones parameters for B are same as Part A). Under Gromacs 4.0.x, you should be able to leave the topology alone entirely. No B state should need to be defined explicitly. I am using default parameters for couple-lambda0 and couple-lambda1. The default behavior is vdw-q, is it not? Then what you're doing still doesn't make sense to me entirely, because no decoupling would even be specified if that is the case. Can you post your actual .mdp file (mdout.mdp might be better if you're leaving it up to grompp to decide the defaults, etc)? I'll admit I haven't done much with the 4.0 free energy code, but it seems to me that with an unaltered topology (no B state explicitly defined) you should be able to set the following in the .mdp file: couple-lambda0 = vdw-q couple-lambda1 = vdw ...to decouple only Coulombic interactions, and likewise: couple-lambda0 = vdw couple-lambda1 = none ...to decouple van der Waals. Perhaps someone can correct me if I've gotten this wrong. For one of my simulation (glucose + ionic liquids) I am getting solvation free energy 35.5 Kcal/mol. Is it too high? Not a clue. Isn't that what you're trying to figure out? Is there any sort of literature value for such a parameter? In any case, it's more plausible than your original value. -Justin Nilesh On Sun, January 10, 2010 1:08 pm, Justin A. Lemkul wrote: On 1/10/10 1:03 PM, Nilesh Dhumal wrote: Hello Justin, Should I add dum_opls no. in atom type file in which Lennard-jones parameters will set to zero. OR. I don't have to modify anything prog.will take care of lennard-jones interactions. I have put charge zero on all atom in TYPE B. You shouldn't have to do make any modifications to the topology whatsoever, according to the documentation. The couple-lambda0 and couple-lambda1 parameters should do the decoupling for you. -Justin Nilesh On Sun, January 3, 2010 12:57 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: THanks Justin. I am using Groamcs 4.0.7 version. I will do more simulation with more lambada value betwen 0 to 0.1. I really think this is less significant than point #1 I mentioned before. Be sure you are actually setting the relevant .mdp parameters (couple-lambda0, couple-lambda1, etc) as they will heavily influence the results. Leaving them as default (as you are doing, per the .mdp file you provided) may give you something other than what you want. And do consult the links I provided; they give excellent background on the nature of these calculations. -Justin Thanks. Nilesh On Sun, January 3, 2010 12:39 pm, Justin A. Lemkul wrote: Nilesh Dhumal wrote: Hey All, I am trying to calculate the solvation free energy of glucose in ionic liquids. I am getting really large value for lamda=0 in vacuum and solvent as well. I am getting final solvation energy value ~ 310 kJ/mol. Can you tell me where I am wrong. There are several possibilities: 1. Are you decoupling LJ and Coulombic interactions simultaneously? They should be done in separate steps, or else you get very large errors. 2. You probably need a lot more lambda points to properly define the curve. For strongly hydrogen-bonding solutes, you should probably include closely-spaced points near lambda = 0 (the fully-coupled state), like 0.005, 0.01, 0.015, etc. The curve from 0 to 0.1 will otherwise be ill-defined and prone to large errors. You also do not say which version of
Re: [gmx-users] Unstable Minimizations
Thanks Justin, I went back to the original pdb files. These were conformations of the same protein derived from molecular dynamics simulations performed by Andrey. What I intially attempted was preping the structures using tleap, hoping to paint in missing atoms for residues. Then use this to replace non-standard residues sed s/PRO\ A\ \ \ 1/NPROA\ \ \ 1/g fzd2_md7-8_c6_cc.pdb | sed s/PRO\ B\ \ \ 1/NPROB\ \ \ 1/g | sed s/PHE\ A\ \ 99/CPHEA\ \ 99/g | sed s/PHE\ B\ \ 99/CPHEB\ \ 99/g | sed s/O\ \ \ CPHE/OC1\ CPHE/g | sed s/OXT\ CPHE/OC2\ CPHE/g | sed s/HIS\ /HID\ /g | sed s/LYS\ /LYP\ /g | sed s/CYS\ /CYN\ /g protein2.pdb Then fixed the nterminal residue name. Finally replaced all CYS to CYS2 I went back and did the same thing except for tleap. It pdb2gmx seems to process these files without needing the tleap step. Still I see some of the same lincs errors. rms 10.669050, max 173.182678 (between atoms 1857 and 1859) rms 10.669803, max 173.177811 (between atoms 1857 and 1859) rms 10.670179, max 173.175400 (between atoms 1857 and 1859) rms 10.670368, max 173.174149 (between atoms 1857 and 1859) rms 10.670460, max 173.173553 (between atoms 1857 and 1859) rms 10.670508, max 173.173141 (between atoms 1857 and 1859) rms 10.670531, max 173.173035 (between atoms 1857 and 1859) rms 10.670543, max 173.172958 (between atoms 1857 and 1859) rms 10.670549, max 173.172928 (between atoms 1857 and 1859) rms 10.670552, max 173.172913 (between atoms 1857 and 1859) rms 10.670554, max 173.172913 (between atoms 1857 and 1859) rms 10.670554, max 173.172913 (between atoms 1857 and 1859) ATOM 1857 CA HIE 120 43.362 28.084 25.727 1.00 0.00 ATOM 1858 HA HIE 120 43.677 27.135 25.748 1.00 0.00 ATOM 1859 CB HIE 120 42.112 28.226 24.788 1.00 0.00 also this atom consistently has a very high Fmax Step=3, Dmax= 1.4e-02 nm, Epot= 1.45860e+10 Fmax= 2.82224e+12, atom= 19392 Step=4, Dmax= 7.2e-03 nm, Epot= 1.45396e+10 Fmax= 2.82207e+12, atom= 19392 Step=5, Dmax= 3.6e-03 nm, Epot= 1.45106e+10 Fmax= 2.82194e+12, atom= 19392 Step=6, Dmax= 1.8e-03 nm, Epot= 1.44953e+10 Fmax= 2.82181e+12, atom= 19392 Step=7, Dmax= 9.0e-04 nm, Epot= 1.44887e+10 Fmax= 2.82196e+12, atom= 19392 Step=8, Dmax= 4.5e-04 nm, Epot= 1.44850e+10 Fmax= 2.82196e+12, atom= 19392 Step=9, Dmax= 2.2e-04 nm, Epot= 1.44832e+10 Fmax= 2.82196e+12, atom= 19392 Step= 10, Dmax= 1.1e-04 nm, Epot= 1.44822e+10 Fmax= 2.82196e+12, atom= 19392 Step= 11, Dmax= 5.6e-05 nm, Epot= 1.44818e+10 Fmax= 2.82196e+12, atom= 19392 Step= 12, Dmax= 2.8e-05 nm, Epot= 1.44815e+10 Fmax= 2.82196e+12, atom= 19392 Step= 13, Dmax= 1.4e-05 nm, Epot= 1.44814e+10 Fmax= 2.82196e+12, atom= 19392 Its not clear to me what we should do to correct this structures...maybe Andrey has some input. http://boinc.drugdiscoveryathome.com/em_restrained_rcs_mdrun2.txt On Sun, Jan 10, 2010 at 5:37 PM, Justin A. Lemkul jalem...@vt.edu wrote: On 1/10/10 5:18 PM, Jack Shultz wrote: I am trying to get this workflow opperational. However, my systems are getting unstable. I have preped two mdp files: 1) one for restrained 2) unrestrained. LINCS errors appear for restrained and unrestrained has infinite energy appearing. http://boinc.drugdiscoveryathome.com/_em_restrained_rcs_mdrun.txt_ http://boinc.drugdiscoveryathome.com/em_restrained_rcs_mdrun.txt __ This log file shows several long bond warnings, which may be the root of your problem. See here: http://www.gromacs.org/Documentation/Errors#Long_bonds_and.2for_missing_atoms Since your minimization is failing immediately, there is something physically unreasonable about your structure, such that EM cannot resolve the problem. Note, too, that one of the long bond warnings pertained to atom 1668, which is the location of the first LINCS warning. Coincidence? Not likely. Re-examine the starting structure and figure out if anything is missing or poorly reconstructed (e.g., from initially missing atoms). This is where I get the LINCS Warnings Step -1, time -0.001 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.461520, max 14.428611 (between atoms 1668 and 1669) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length Steepest Descents: Tolerance (Fmax) = 1.0e+04 Number of steps= 100 Warning: 1-4 interaction between 1658 and 1672 at distance 2.655 which is larger than the 1-4 table size 2.400 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size Step=0, Dmax= 1.0e-02 nm, Epot= 1.09364e+09 Fmax= 2.21154e+11, atom= 3292 Step 1, time 0.001 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.680509, max 22.293625 (between atoms 1668 and 1670) bonds that rotated more than 30 degrees: atom 1 atom 2
Re: [gmx-users] Unstable Minimizations
On 1/10/10 9:47 PM, Jack Shultz wrote: Thanks Justin, I went back to the original pdb files. These were conformations of the same protein derived from molecular dynamics simulations performed by Andrey. What I intially attempted was preping the structures using tleap, hoping to paint in missing atoms for residues. Then use this to replace Well, it seems that you may be hoping for too much :) Your log file shows a whole bunch of failures that look to be related to some early processing of your structure, and other warnings about close contacts detected in tleap. I think you may need to start with an actual intact structure, or else coax your preparation steps to make this happen. I am not too familiar with tleap and sleap, do they magically fix missing atoms? -Justin non-standard residues sed s/PRO\ A\ \ \ 1/NPROA\ \ \ 1/g fzd2_md7-8_c6_cc.pdb | sed s/PRO\ B\ \ \ 1/NPROB\ \ \ 1/g | sed s/PHE\ A\ \ 99/CPHEA\ \ 99/g | sed s/PHE\ B\ \ 99/CPHEB\ \ 99/g | sed s/O\ \ \ CPHE/OC1\ CPHE/g | sed s/OXT\ CPHE/OC2\ CPHE/g | sed s/HIS\ /HID\ /g | sed s/LYS\ /LYP\ /g | sed s/CYS\ /CYN\ /g protein2.pdb Then fixed the nterminal residue name. Finally replaced all CYS to CYS2 I went back and did the same thing except for tleap. It pdb2gmx seems to process these files without needing the tleap step. Still I see some of the same lincs errors. rms 10.669050, max 173.182678 (between atoms 1857 and 1859) rms 10.669803, max 173.177811 (between atoms 1857 and 1859) rms 10.670179, max 173.175400 (between atoms 1857 and 1859) rms 10.670368, max 173.174149 (between atoms 1857 and 1859) rms 10.670460, max 173.173553 (between atoms 1857 and 1859) rms 10.670508, max 173.173141 (between atoms 1857 and 1859) rms 10.670531, max 173.173035 (between atoms 1857 and 1859) rms 10.670543, max 173.172958 (between atoms 1857 and 1859) rms 10.670549, max 173.172928 (between atoms 1857 and 1859) rms 10.670552, max 173.172913 (between atoms 1857 and 1859) rms 10.670554, max 173.172913 (between atoms 1857 and 1859) rms 10.670554, max 173.172913 (between atoms 1857 and 1859) ATOM 1857 CA HIE 120 43.362 28.084 25.727 1.00 0.00 ATOM 1858 HA HIE 120 43.677 27.135 25.748 1.00 0.00 ATOM 1859 CB HIE 120 42.112 28.226 24.788 1.00 0.00 also this atom consistently has a very high Fmax Step=3, Dmax= 1.4e-02 nm, Epot= 1.45860e+10 Fmax= 2.82224e+12, atom= 19392 Step=4, Dmax= 7.2e-03 nm, Epot= 1.45396e+10 Fmax= 2.82207e+12, atom= 19392 Step=5, Dmax= 3.6e-03 nm, Epot= 1.45106e+10 Fmax= 2.82194e+12, atom= 19392 Step=6, Dmax= 1.8e-03 nm, Epot= 1.44953e+10 Fmax= 2.82181e+12, atom= 19392 Step=7, Dmax= 9.0e-04 nm, Epot= 1.44887e+10 Fmax= 2.82196e+12, atom= 19392 Step=8, Dmax= 4.5e-04 nm, Epot= 1.44850e+10 Fmax= 2.82196e+12, atom= 19392 Step=9, Dmax= 2.2e-04 nm, Epot= 1.44832e+10 Fmax= 2.82196e+12, atom= 19392 Step= 10, Dmax= 1.1e-04 nm, Epot= 1.44822e+10 Fmax= 2.82196e+12, atom= 19392 Step= 11, Dmax= 5.6e-05 nm, Epot= 1.44818e+10 Fmax= 2.82196e+12, atom= 19392 Step= 12, Dmax= 2.8e-05 nm, Epot= 1.44815e+10 Fmax= 2.82196e+12, atom= 19392 Step= 13, Dmax= 1.4e-05 nm, Epot= 1.44814e+10 Fmax= 2.82196e+12, atom= 19392 Its not clear to me what we should do to correct this structures...maybe Andrey has some input. http://boinc.drugdiscoveryathome.com/em_restrained_rcs_mdrun2.txt On Sun, Jan 10, 2010 at 5:37 PM, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: On 1/10/10 5:18 PM, Jack Shultz wrote: I am trying to get this workflow opperational. However, my systems are getting unstable. I have preped two mdp files: 1) one for restrained 2) unrestrained. LINCS errors appear for restrained and unrestrained has infinite energy appearing. http://boinc.drugdiscoveryathome.com/_em_restrained_rcs_mdrun.txt_ http://boinc.drugdiscoveryathome.com/em_restrained_rcs_mdrun.txt __ This log file shows several long bond warnings, which may be the root of your problem. See here: http://www.gromacs.org/Documentation/Errors#Long_bonds_and.2for_missing_atoms Since your minimization is failing immediately, there is something physically unreasonable about your structure, such that EM cannot resolve the problem. Note, too, that one of the long bond warnings pertained to atom 1668, which is the location of the first LINCS warning. Coincidence? Not likely. Re-examine the starting structure and figure out if anything is missing or poorly reconstructed (e.g., from initially missing atoms). This is where I get the LINCS Warnings Step -1, time -0.001 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.461520, max 14.428611 (between atoms 1668 and 1669) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length Steepest Descents:
