Re: [gmx-users] keeping water (entirely) out of the bilayer core

2013-05-05 Thread XAvier Periole

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)

2013-05-05 Thread Group Gro
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)

2013-05-05 Thread Justin Lemkul



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

2013-05-05 Thread Kong xq
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

2013-05-05 Thread Shima Arasteh
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

2013-05-05 Thread Justin Lemkul



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

2013-05-05 Thread Martin Siegert
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

2013-05-05 Thread Christopher Neale
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

2013-05-05 Thread Christopher Neale
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

2013-05-05 Thread Mark Abraham
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

2013-05-05 Thread Mark Abraham
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

2013-05-05 Thread Jesper Sørensen
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

2013-05-05 Thread Kong xq
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

2013-05-05 Thread shahid nayeem
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