Re: [gmx-users] keeping water (entirely) out of the bilayer core
You could define a repulsive potential that would apply to water molecules. I am not sure how it is called but it is available and described in the manual. On May 4, 2013, at 2:46, Christopher Neale chris.ne...@mail.utoronto.ca wrote: Dear users: I am interested in running simulations of lipid bilayers in which I keep all water molecules out of the bilayer core (not just statistically, but absolutely). However, I have been unable to figure out how to do it. I'll list what I have tried in the hope that others have some ideas or even perhaps know how to do this with standard gromacs. Everything that I have tried revolves around using the pull code, setting the entire lipid bilayer as the reference group, and having thousands of pulled groups -- one for each water molecule. Also, I modified the gromacs source code in mdlib/pull.c (version 4.6.1) so that the force is only applied when the displacement is smaller than the desired distance and not when the displacement is larger than the specified distance (to keep water out but then to otherwise let it go anywhere without bias): static void do_pull_pot(int ePull, t_pull *pull, t_pbc *pbc, double t, real lambda, real *V, tensor vir, real *dVdl) { intg, j, m; dvec dev; double ndr, invdr; real k, dkdl; t_pullgrp *pgrp; /* loop over the groups that are being pulled */ *V= 0; *dVdl = 0; for (g = 1; g 1+pull-ngrp; g++) { pgrp = pull-grp[g]; get_pullgrp_distance(pull, pbc, g, t, pgrp-dr, dev); /* ## HERE IS THE CODE CHANGE if(dev[0]0.0){ dev[0]=0.0; } /* End code snippit */ The most relevant lines from the .mdp file were: pull = umbrella pull_geometry= distance pull_dim = N N Y pull_start = no pull_ngroups = 9000 pull_group0 = POPC pull_group1 = r_1__OW pull_init1 = 1.9 pull_k1 = 500 pull_group2 = r_2__OW pull_init2 = 1.9 pull_k2 = 500 ... etc... The problem with this was that it crashed immediately with an error that my pulling distance was greater than 1/2 of the box vector.: Fatal error: Distance of pull group 165 (5.353939 nm) is larger than 0.49 times the box size (5.395960) Pretty obvious in retrospect. If I could get this error to go away, then everything should be fine because I have modified the code so that forces are zero at large displacements. However, I am worried that if I simply modify grompp to bypass this error then I might get some type of strange and possibly silent error in mdrun. Any thoughts on this are really appreciated. For my second try, I used pull_geometry = direction-periodic (see more details below), but that also didn't work because if I set pull_vec1 = 0 0 1, then everything gets pulled to the upper part of the box (and LINCS errors ensue), rather than simply pulling away from the bilayer. Pcoupl= berendsen pcoupltype = semiisotropic compressibility= 4.5e-5 0 pull = umbrella pull_geometry= direction-periodic pull_start = no pull_ngroups = 9000 pull_group0 = POPC pull_group1 = r_1__OW pull_init1 = 1.9 pull_k1 = 500 pull_vec1= 0 0 1 pull_group2 = r_2__OW pull_init2 = 1.9 pull_k2 = 500 pull_vec2= 0 0 1 ... etc... If anybody has an idea about how I could get this to work with standard gromacs or with a modified version, I would be very grateful to hear your thoughts. Thank you, 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/Support/Mailing_Lists/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/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists
[gmx-users] (no subject)
Dear GROMACS users, I am working on protein-ligand complexes and when I run mdrun -deffnm nvt.tpr -v I run into this error: Fatal error: 2 particles communicated to PME node 1 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x. I found the best docked position of my ligand by Autodock run and copied the best position of my ligand in a pdb format. Then I ran PRODRG to provide gro and itp files. I have tried different drugs as ligand and all of them are OK. I searched the mailing list and found other users have had this error. I understood that my system would be unstable. What should I do to solve this problem? Thanks alot. S. Faraji -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists
Re: [gmx-users] (no subject)
On 5/5/13 5:10 AM, Group Gro wrote: Dear GROMACS users, I am working on protein-ligand complexes and when I run mdrun -deffnm nvt.tpr -v I run into this error: Fatal error: 2 particles communicated to PME node 1 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x. I found the best docked position of my ligand by Autodock run and copied the best position of my ligand in a pdb format. Then I ran PRODRG to provide gro and itp files. I have tried different drugs as ligand and all of them are OK. I searched the mailing list and found other users have had this error. I understood that my system would be unstable. What should I do to solve this problem? Start with a better topology. PRODRG topologies produce bad results. http://pubs.acs.org/doi/abs/10.1021/ci100335w -Justin -- Justin A. Lemkul, Ph.D. Research Scientist 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/Support/Mailing_Lists/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/Support/Mailing_Lists
[gmx-users] REMD Statistics
Dear GMX users, I have some concerns about the statistics analysis of REMD which do need your generous help. I performed a 50ns isothermal-isobaric REMD simulation with 64 replicas spaning from 300K to 390K. Then I want to do some statistics analysis for the results. First, I calculated the acceptance ratio for the adjacent pairs of temperatures (showing below). It seems that most of the values are uniform ( 0.25 +/- 0.11) . Second, I followed the track of replica 1 in the temperature space with the information from replica_temp.xvg generated by demux.pl(demux.pl md0.log) scipt. However, the histogram(table below) indicated that the distribution of replica 1 in temperature space is not very uniform. So my question is why did the nonuniform temperature distribution of replica 1 occur given the fact that the acceptance ratio indicated a free random walk in temperature space ? What's more, whether a longer simulation was needed to make a more flatter temperature distribution for each replica? By the way, I want to double check that each column except the first in replica_temp.xvg represents a constant temperature corresponding to that specified in tpr file, we may need to follow the replica index (such as 0 for the first replica) among different coloumns (or temperatures). Does this make sense? Best regards! Xiangqian Kong *The frequency of replica 1 in temperature space * TemperatureFrequency31016.91%32011.58%33013.68%34013.79%35012.64%36011.60% 3704.38%3805.26%39010.16% -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists
[gmx-users] RMSD
Hi, I' like to know if it is possible to get the average RMSD through g_rms command? Or I need to get it manually? Thanks for your suggestions. Sincerely, Shima -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists
Re: [gmx-users] RMSD
On 5/5/13 12:40 PM, Shima Arasteh wrote: Hi, I' like to know if it is possible to get the average RMSD through g_rms command? Or I need to get it manually? Use g_analyze. -Justin -- Justin A. Lemkul, Ph.D. Research Scientist 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/Support/Mailing_Lists/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/Support/Mailing_Lists
Re: [gmx-users] configuration/installation of gromacs-4.6.1 on heterogeneous cluster
Hi, On Sat, May 04, 2013 at 09:50:37AM +0200, Mark Abraham wrote: On Sat, May 4, 2013 at 1:41 AM, Martin Siegert sieg...@sfu.ca wrote: Hi, I am struggling with the configuration and compilation/installation of gromacs-4.6.1. Our cluster has 2 different processors: the older generation supports sse4.1, the newer sse4.2. Configuration and compilation must be done on the headnode of the cluster, which supports sse4.2. I am using the following command to configure gromacs-4.6.1: CFLAGS='-fpic -O3 -axSSE4.2,SSE4.1 -xSSSE3 -ip -opt-prefetch' \ CXXFLAGS='-fpic -O3 -axSSE4.2,SSE4.1 -xSSSE3 -ip -opt-prefetch' \ FFLAGS='-fpic -O3 -axSSE4.2,SSE4.1 -xSSSE3 -ip -opt-prefetch' \ If the above is successful at generating instructions for SSE4.1... -xSSSE3 ensures that the resulting code will run on processors that support ssse3 or later and is basically equivalent to -mssse3. -axSSE4.2,SSE4.1 tells the compiler to generate multiple code paths for sse4.1 and sse4.2 while still supporting the default architecture ssse3 i.e., when run on a processor that supports sse4.2 then the code will use the special instructions for that architecture). We compile all our codes with -axSSE4.2,SSE4.1 -xSSSE3 and they run on ssse3 and better. CC=mpicc \ CXX=mpicxx \ FC=mpif90 \ LDFLAGS=-lfftw3f -lgoto2 -Wl,-rpath,/usr/local/gromacs-4.6.1/lib \ cmake -DCMAKE_VERBOSE_MAKEFILE=ON -DGMX_MPI=ON \ -DGMX_CPU_ACCELERATION=SSE4.1 -DGMX_OPENMP=OFF -DGMX_GPU=OFF \ -DCMAKE_INSTALL_PREFIX=/usr/local/gromacs-4.6.1 \ -DCMAKE_SKIP_RPATH=YES \ ../gromacs-4.6.1 However, after compilation/installation mdrun_mpi fails on nodes that only support sse4.1 with Illegal instruction. ... then this is not a surprise. Your compiler has been allowed to generate SSE 4.2 instructions, I suspect. Not really ... -xssse3 specifies the default architecture that the code supports. The CMakeCache.txt file contains the line: BUILD_CPU_FEATURES:INTERNAL=aes apic clfsh cmov cx8 cx16 htt lahf_lm mmx msr nonstop_tsc pcid pclmuldq pdcm pdpe1gb popcnt pse rdtscp sse2 sse3 sse4.1 sse4.2 ssse3 Since this line contains sse4.2 it appears that the flag -DGMX_CPU_ACCELERATION=SSE4.1 is ignored. That list is what features are available on the CPU (mostly for helping detect what acceleration to use, and to solve problems). The target acceleration scheme (i.e. code with heavy use of compiler instrinsics) is SSE4.1, which is what GROMACS will compile for if you leave it alone :-) SSE4.2 is roughly useless for mdrun. What is the correct way of specifying the cpu architecture within the cmake build system? (I never had problems with this with the pre 4.6 versions). Back then, GROMACS was nearly insensitive to compilation settings, because of the use of assembly kernels. Now GROMACS is sensitive to compiler version (in that compilation of SIMD instrinsics needs to work well, and OpenMP needs to be supported, etc.) but one still doesn't generally want to mess with the compiler flags. We have some internal disgreement about whether we should be permitting/encouraging/facilitiating setting/checking compiler flags. So far nobody has demonstrated a use case that suggests we need to support more than shut up and let GROMACS do its thing! For curiosity I tried this, i.e., I did not set CFLAGS, CXXFLAGS and FFLAGS and use the gromacs supplied values which turn out to be -O3 -ip -fPIC -msse4.1 The result is absolutely the same: on the nodes that only support sse4.1 the code fails with Illegal instruction. Any suggestions what else I should try? I am mostly concerned about the following setting in CMakeCache.txt: //Build CPU brand BUILD_CPU_BRAND:INTERNAL=Intel(R) Xeon(R) CPU X5650 @ 2.67GHz //Build CPU family BUILD_CPU_FAMILY:INTERNAL=6 //Build CPU features BUILD_CPU_FEATURES:INTERNAL=aes apic clfsh cmov cx8 cx16 htt lahf_lm mmx msr nonstop_tsc pcid pclmuldq pdcm pdpe1gb popcnt pse rdtscp sse2 sse3 sse4.1 sse4.2 ssse3 //Build CPU model BUILD_CPU_MODEL:INTERNAL=44 //Build CPU stepping BUILD_CPU_STEPPING:INTERNAL=2 As the names of the variables say, these are for the BUILD cpu, not necessarily for the cpu the code needs to run on. Thus, if any of these are used to generate processor specific code, then this would explain the illegal instruction. Consequently I would need to change these settings, but to what? Cheers, Martin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists
[gmx-users] keeping water (entirely) out of the bilayer core
Thank you XAvier, The only thing I could find in the 4.6.1 manual about repulsive interactions is the repulsive part of the LJ and Buckingham potentials, in addition to tabulated potentials. Were you referring to one of these? Chris. -- original message -- You could define a repulsive potential that would apply to water molecules. I am not sure how it is called but it is available and described in the manual. On May 4, 2013, at 2:46, Christopher Neale chris.neale at mail.utoronto.ca wrote: Dear users: I am interested in running simulations of lipid bilayers in which I keep all water molecules out of the bilayer core (not just statistically, but absolutely). However, I have been unable to figure out how to do it. I'll list what I have tried in the hope that others have some ideas or even perhaps know how to do this with standard gromacs. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists
[gmx-users] keeping water (entirely) out of the bilayer core
Thank you for the advice Jianguo. This seems like a good way forward. I was hoping not to have to learn how to use new software, but perhaps it is time. Thank you, Chris. -- original message -- Hi Chris, Perhaps you can try PLUMED+GROMACS. In that case, you can define a collective variable as mindist between the water molecules and the lipid tails. And then apply wall potential to keep this CV above a certain value. Jianguo Dear users: I am interested in running simulations of lipid bilayers in which I keep all water molecules out of the bilayer core(not just statistically, but absolutely). However, I have been unable to figure out how to do it. I'll list what I have tried in the hope that others have some ideas or even perhaps know how to do this with standard gromacs. ... -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists
Re: [gmx-users] configuration/installation of gromacs-4.6.1 on heterogeneous cluster
On Sun, May 5, 2013 at 8:12 PM, Martin Siegert sieg...@sfu.ca wrote: Hi, On Sat, May 04, 2013 at 09:50:37AM +0200, Mark Abraham wrote: On Sat, May 4, 2013 at 1:41 AM, Martin Siegert sieg...@sfu.ca wrote: Hi, I am struggling with the configuration and compilation/installation of gromacs-4.6.1. Our cluster has 2 different processors: the older generation supports sse4.1, the newer sse4.2. Configuration and compilation must be done on the headnode of the cluster, which supports sse4.2. I am using the following command to configure gromacs-4.6.1: CFLAGS='-fpic -O3 -axSSE4.2,SSE4.1 -xSSSE3 -ip -opt-prefetch' \ CXXFLAGS='-fpic -O3 -axSSE4.2,SSE4.1 -xSSSE3 -ip -opt-prefetch' \ FFLAGS='-fpic -O3 -axSSE4.2,SSE4.1 -xSSSE3 -ip -opt-prefetch' \ If the above is successful at generating instructions for SSE4.1... -xSSSE3 ensures that the resulting code will run on processors that support ssse3 or later and is basically equivalent to -mssse3. -axSSE4.2,SSE4.1 tells the compiler to generate multiple code paths for sse4.1 and sse4.2 while still supporting the default architecture ssse3 i.e., when run on a processor that supports sse4.2 then the code will use the special instructions for that architecture). OK - but still not something from which GROMACS can derive any significant advantage, even if it works reliably :-) We compile all our codes with -axSSE4.2,SSE4.1 -xSSSE3 and they run on ssse3 and better. CC=mpicc \ CXX=mpicxx \ FC=mpif90 \ LDFLAGS=-lfftw3f -lgoto2 -Wl,-rpath,/usr/local/gromacs-4.6.1/lib \ cmake -DCMAKE_VERBOSE_MAKEFILE=ON -DGMX_MPI=ON \ -DGMX_CPU_ACCELERATION=SSE4.1 -DGMX_OPENMP=OFF -DGMX_GPU=OFF \ -DCMAKE_INSTALL_PREFIX=/usr/local/gromacs-4.6.1 \ -DCMAKE_SKIP_RPATH=YES \ ../gromacs-4.6.1 However, after compilation/installation mdrun_mpi fails on nodes that only support sse4.1 with Illegal instruction. ... then this is not a surprise. Your compiler has been allowed to generate SSE 4.2 instructions, I suspect. Not really ... -xssse3 specifies the default architecture that the code supports. The CMakeCache.txt file contains the line: BUILD_CPU_FEATURES:INTERNAL=aes apic clfsh cmov cx8 cx16 htt lahf_lm mmx msr nonstop_tsc pcid pclmuldq pdcm pdpe1gb popcnt pse rdtscp sse2 sse3 sse4.1 sse4.2 ssse3 Since this line contains sse4.2 it appears that the flag -DGMX_CPU_ACCELERATION=SSE4.1 is ignored. That list is what features are available on the CPU (mostly for helping detect what acceleration to use, and to solve problems). The target acceleration scheme (i.e. code with heavy use of compiler instrinsics) is SSE4.1, which is what GROMACS will compile for if you leave it alone :-) SSE4.2 is roughly useless for mdrun. What is the correct way of specifying the cpu architecture within the cmake build system? To clarify: setting GMX_CPU_ACCELERATION=SSE4.1 is all that you require. IIRC the build-cpu detection still runs in this case, but the results have no real bearing on the outcome. You can see the value of that variable in the CMakeCache.txt. The -msse4.1 compiler flag is used if and only if that is the value of that variable. (I never had problems with this with the pre 4.6 versions). Back then, GROMACS was nearly insensitive to compilation settings, because of the use of assembly kernels. Now GROMACS is sensitive to compiler version (in that compilation of SIMD instrinsics needs to work well, and OpenMP needs to be supported, etc.) but one still doesn't generally want to mess with the compiler flags. We have some internal disgreement about whether we should be permitting/encouraging/facilitiating setting/checking compiler flags. So far nobody has demonstrated a use case that suggests we need to support more than shut up and let GROMACS do its thing! For curiosity I tried this, i.e., I did not set CFLAGS, CXXFLAGS and FFLAGS and use the gromacs supplied values which turn out to be -O3 -ip -fPIC -msse4.1 The result is absolutely the same: on the nodes that only support sse4.1 the code fails with Illegal instruction. That is extremely surprising. As above, the appearance of -msse4.1 indicates GROMACS is doing what is expected for your case. Any suggestions what else I should try? Compiling on the sse4.1 node would be an interesting experiment, but you seem to suggest this is impossible. Otherwise, this could be a compiler bug. GROMACS has a long history of pushing compilers too far ;-) For example, we know intel 11.1 has some bug with SSE4.1 code generation ( http://redmine.gromacs.org/issues/1126), but whatever is the cause here is a different issue. What is the full compiler version? I am mostly concerned about the following setting in CMakeCache.txt: //Build CPU brand BUILD_CPU_BRAND:INTERNAL=Intel(R) Xeon(R) CPU X5650 @ 2.67GHz //Build CPU family
Re: [gmx-users] REMD Statistics
On Sun, May 5, 2013 at 5:14 PM, Kong xq xqkong...@gmail.com wrote: Dear GMX users, I have some concerns about the statistics analysis of REMD which do need your generous help. I performed a 50ns isothermal-isobaric REMD simulation with 64 replicas spaning from 300K to 390K. Then I want to do some statistics analysis for the results. First, I calculated the acceptance ratio for the adjacent pairs of temperatures (showing below). It seems that most of the values are uniform ( 0.25 +/- 0.11) . That is an abnormally high variation in my experience. That could suggest lack of generalized equilibrium on this time scale. Second, I followed the track of replica 1 in the temperature space with the information from replica_temp.xvg generated by demux.pl(demux.pl md0.log) scipt. However, the histogram(table below) indicated that the distribution of replica 1 in temperature space is not very uniform. So my question is why did the nonuniform temperature distribution of replica 1 occur given the fact that the acceptance ratio indicated a free random walk in temperature space ? The acceptance ratio is only a proxy for the existence of a random walk. It would be an amusing exercise to construct a toy system with replicas that would only exchange over a few adjacent temperatures. With suitable overlap of replicas constrained-temperature subranges, you could synthesise uniform exchange rates with no interesting replica flow. For example, if * state A is possible in the lower two-thirds of your temperature range, and * state B is possible in the upper two-thirds of your temperature range, and * no other states are possible, and * your initial configurations include examples of both A and B, and * A and B have no available interconversion pathways, then you can see apparently uniform neighbour exchange rates, and the REMD simulation is probably not being very useful in the normal case of trying to enhance statistics at the lowest temperature. The possibility of this kind of kinetic trapping is discussed in the REMD literature (including some rather old mentions). The extreme challenge of actually achieving replica flow efficiently is discussed in a series of papers from Nadler and Hansmann. Even spacing the temperatures for even replica flow for met-enkephalin is not trivial. What's more, whether a longer simulation was needed to make a more flatter temperature distribution for each replica? Probably, but as above, not necessarily meaningful. It might be instructive to construct some random walks following the GROMACS REMD protocal with acceptance rate 0.25 and see whether your intuition is valid :-) By the way, I want to double check that each column except the first in replica_temp.xvg represents a constant temperature corresponding to that specified in tpr file, we may need to follow the replica index (such as 0 for the first replica) among different coloumns (or temperatures). Does this make sense? Look at the first row and column, and back at the .log files to see what exchanges take place. That's much better than trying a written description! Mark Best regards! Xiangqian Kong *The frequency of replica 1 in temperature space * TemperatureFrequency31016.91%32011.58%33013.68%34013.79%35012.64%36011.60% 3704.38%3805.26%39010.16% -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists
Re: [gmx-users] constant protonation state MD
Shahid, This must be system dependent then because there are experimental methods, e.g. NMR and CD, that can help determine the secondary structure of proteins/peptides are various pH ranges. What is your system? A peptide? A globular protein? Best, Jesper On May 3, 2013, at 7:33 PM, shahid nayeem msnay...@gmail.com wrote: Thanks a lot Erik and Baptista I am interested in simulating the change in secondary structure which is supposed to be influenced by the change in the pH environment of the cell. Experimentally it is not known but it has been proposed by many that the change in pH leads to change in conformation. I did constant protonation state MD at two different pH. I got an alpha helix converting to beta sheet at higher pH but not at lower pH. I am bothered to prove the reliability of the simulation results as experimentally it cant be established. Shahid On Sat, May 4, 2013 at 5:46 AM, bapti...@itqb.unl.pt wrote: Hi Shahid, As Erik said, it depends... on your system, on the process you are studying in that system, on the property you think it's relevant to study that process, etc. If your question refers to the (de)protonation of acidic and basic groups usually occuring in aqueous solution, there are methods to perform what is called constant-pH MD, where the protonation states of those groups change in a discrete or continuous way along the simulation. If you are interested, start with the corresponding GROMACS how-to ( http://www.gromacs.org/**Documentation/How-tos/**Constant_pH_Simulationhttp://www.gromacs.org/Documentation/How-tos/Constant_pH_Simulation) and then search the literature. In any case, you don't need any of that fancy stuff unless you have reasons to think that the properties you want to study in your system are strongly dependent on protonation-conformation coupling effects. In some cases you may be suspicious about the effect of the protonation state of one group (e.g., a histidine), but that can be simply tested by performing two constant-protonation MD simulations, one for each state (you can even try to use a linear response approximation on top of that). But in most cases you don't need even that. For what it's worth: I was one of the first to develop and apply constant-pH MD and I don't bother to use it most of the time... :) Best, Antonio On Fri, 3 May 2013, Erik Marklund wrote: Yes that's what lambda dynamics does. I mentioned it since it addresses the interplay between protonation and structure. So to answer your original question: it depends. Erik On 3 May 2013, at 15:27, shahid nayeem msnay...@gmail.com wrote: If I know correctly in lambda dynamics the dynamics of protonation/deprotonation equilibria is accounted for while my question relates to the typical constant protonation MD where each titratable group remains in one protonation state throughout the simulation. Please educate me Shahid On Fri, May 3, 2013 at 6:22 PM, Erik Marklund er...@xray.bmc.uu.se wrote: I don't have one in mind. It's a delicate question and perhaps I shouldn't have phrased it the way I did. Nevertheless, the pKa of most side chains mean that their protonation will be dominated by one state for most pH values. pKa-shifts and complicated interplay between protonation and structure cause exceptions to this and you should be aware that such things may in some cases be important. Also consider the timescales of pH-depedent structural changes and the length of your simulations. You could check out the papers on lambda dynamics by C. Brooks III for an interesting take on sampling multiple protonation states. Best, Erik On 3 May 2013, at 14:05, shahid nayeem msnay...@gmail.com wrote: Thanks a lot Erik. Could I get some reference based on which you say that much of the structural biology will be largely unaffected. Shahid On Fri, May 3, 2013 at 1:05 PM, Erik Marklund er...@xray.bmc.uu.se wrote: There's no general answer to that. Proton conductivity measurements, for instance, will be horribly wrong without dynamic protonation. Much (but not all) structural biology, however, will be largely unaffected. Erik On 3 May 2013, at 04:30, shahid nayeem msnay...@gmail.com wrote: Dear all Can someone enlighten me on the reliability of the results obtained from constant protonation state (assigned by different pKa value at different pH) MD simulation. Also want to know its reliability in case of implicit solvation model such as PB/GB calculation. Shahid -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/**mailman/listinfo/gmx-usershttp://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/**Support/Mailing_Lists/Searchhttp://www.gromacs.org/Support/Mailing_Lists/Searchbefore posting! * Please don't post (un)subscribe requests to the list. Use the www
Re: [gmx-users] REMD Statistics
Hi Mark, Thanks for your great help. I am sorry for the negligence to state the variation value correctly( it should be 0.011 rather than 0.11). Does this somewhat small value indicate the generalized equilibrium achieved? I will search the papers you suggested. I am wondering whether the histogram is the gold standard for effectively constructing REMD flow. Moreover, how to do the potential energy overlap analysis in Gromacs, from the edr file for each replica ?Does any tool avalible in Gromacs to do this job ? Xiangqian Kong 2013/5/6 Mark Abraham mark.j.abra...@gmail.com On Sun, May 5, 2013 at 5:14 PM, Kong xq xqkong...@gmail.com wrote: Dear GMX users, I have some concerns about the statistics analysis of REMD which do need your generous help. I performed a 50ns isothermal-isobaric REMD simulation with 64 replicas spaning from 300K to 390K. Then I want to do some statistics analysis for the results. First, I calculated the acceptance ratio for the adjacent pairs of temperatures (showing below). It seems that most of the values are uniform ( 0.25 +/- 0.11) . That is an abnormally high variation in my experience. That could suggest lack of generalized equilibrium on this time scale. Second, I followed the track of replica 1 in the temperature space with the information from replica_temp.xvg generated by demux.pl(demux.pl md0.log) scipt. However, the histogram(table below) indicated that the distribution of replica 1 in temperature space is not very uniform. So my question is why did the nonuniform temperature distribution of replica 1 occur given the fact that the acceptance ratio indicated a free random walk in temperature space ? The acceptance ratio is only a proxy for the existence of a random walk. It would be an amusing exercise to construct a toy system with replicas that would only exchange over a few adjacent temperatures. With suitable overlap of replicas constrained-temperature subranges, you could synthesise uniform exchange rates with no interesting replica flow. For example, if * state A is possible in the lower two-thirds of your temperature range, and * state B is possible in the upper two-thirds of your temperature range, and * no other states are possible, and * your initial configurations include examples of both A and B, and * A and B have no available interconversion pathways, then you can see apparently uniform neighbour exchange rates, and the REMD simulation is probably not being very useful in the normal case of trying to enhance statistics at the lowest temperature. The possibility of this kind of kinetic trapping is discussed in the REMD literature (including some rather old mentions). The extreme challenge of actually achieving replica flow efficiently is discussed in a series of papers from Nadler and Hansmann. Even spacing the temperatures for even replica flow for met-enkephalin is not trivial. What's more, whether a longer simulation was needed to make a more flatter temperature distribution for each replica? Probably, but as above, not necessarily meaningful. It might be instructive to construct some random walks following the GROMACS REMD protocal with acceptance rate 0.25 and see whether your intuition is valid :-) By the way, I want to double check that each column except the first in replica_temp.xvg represents a constant temperature corresponding to that specified in tpr file, we may need to follow the replica index (such as 0 for the first replica) among different coloumns (or temperatures). Does this make sense? Look at the first row and column, and back at the .log files to see what exchanges take place. That's much better than trying a written description! Mark Best regards! Xiangqian Kong *The frequency of replica 1 in temperature space * TemperatureFrequency31016.91%32011.58%33013.68%34013.79%35012.64%36011.60% 3704.38%3805.26%39010.16% -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/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/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! *
Re: [gmx-users] constant protonation state MD
It is a peptide. Shahid On Mon, May 6, 2013 at 4:46 AM, Jesper Sørensen jesoren...@ucsd.edu wrote: Shahid, This must be system dependent then because there are experimental methods, e.g. NMR and CD, that can help determine the secondary structure of proteins/peptides are various pH ranges. What is your system? A peptide? A globular protein? Best, Jesper On May 3, 2013, at 7:33 PM, shahid nayeem msnay...@gmail.com wrote: Thanks a lot Erik and Baptista I am interested in simulating the change in secondary structure which is supposed to be influenced by the change in the pH environment of the cell. Experimentally it is not known but it has been proposed by many that the change in pH leads to change in conformation. I did constant protonation state MD at two different pH. I got an alpha helix converting to beta sheet at higher pH but not at lower pH. I am bothered to prove the reliability of the simulation results as experimentally it cant be established. Shahid On Sat, May 4, 2013 at 5:46 AM, bapti...@itqb.unl.pt wrote: Hi Shahid, As Erik said, it depends... on your system, on the process you are studying in that system, on the property you think it's relevant to study that process, etc. If your question refers to the (de)protonation of acidic and basic groups usually occuring in aqueous solution, there are methods to perform what is called constant-pH MD, where the protonation states of those groups change in a discrete or continuous way along the simulation. If you are interested, start with the corresponding GROMACS how-to ( http://www.gromacs.org/**Documentation/How-tos/**Constant_pH_Simulation http://www.gromacs.org/Documentation/How-tos/Constant_pH_Simulation) and then search the literature. In any case, you don't need any of that fancy stuff unless you have reasons to think that the properties you want to study in your system are strongly dependent on protonation-conformation coupling effects. In some cases you may be suspicious about the effect of the protonation state of one group (e.g., a histidine), but that can be simply tested by performing two constant-protonation MD simulations, one for each state (you can even try to use a linear response approximation on top of that). But in most cases you don't need even that. For what it's worth: I was one of the first to develop and apply constant-pH MD and I don't bother to use it most of the time... :) Best, Antonio On Fri, 3 May 2013, Erik Marklund wrote: Yes that's what lambda dynamics does. I mentioned it since it addresses the interplay between protonation and structure. So to answer your original question: it depends. Erik On 3 May 2013, at 15:27, shahid nayeem msnay...@gmail.com wrote: If I know correctly in lambda dynamics the dynamics of protonation/deprotonation equilibria is accounted for while my question relates to the typical constant protonation MD where each titratable group remains in one protonation state throughout the simulation. Please educate me Shahid On Fri, May 3, 2013 at 6:22 PM, Erik Marklund er...@xray.bmc.uu.se wrote: I don't have one in mind. It's a delicate question and perhaps I shouldn't have phrased it the way I did. Nevertheless, the pKa of most side chains mean that their protonation will be dominated by one state for most pH values. pKa-shifts and complicated interplay between protonation and structure cause exceptions to this and you should be aware that such things may in some cases be important. Also consider the timescales of pH-depedent structural changes and the length of your simulations. You could check out the papers on lambda dynamics by C. Brooks III for an interesting take on sampling multiple protonation states. Best, Erik On 3 May 2013, at 14:05, shahid nayeem msnay...@gmail.com wrote: Thanks a lot Erik. Could I get some reference based on which you say that much of the structural biology will be largely unaffected. Shahid On Fri, May 3, 2013 at 1:05 PM, Erik Marklund er...@xray.bmc.uu.se wrote: There's no general answer to that. Proton conductivity measurements, for instance, will be horribly wrong without dynamic protonation. Much (but not all) structural biology, however, will be largely unaffected. Erik On 3 May 2013, at 04:30, shahid nayeem msnay...@gmail.com wrote: Dear all Can someone enlighten me on the reliability of the results obtained from constant protonation state (assigned by different pKa value at different pH) MD simulation. Also want to know its reliability in case of implicit solvation model such as PB/GB calculation. Shahid -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/**mailman/listinfo/gmx-users http://lists.gromacs.org/mailman/listinfo/gmx-users * Please