Hmm. I am getting the desired output now, at least as far as the pulled atoms are concerned, they are all close to the bilayer center (was not so in a free simulation), and therefore my suspicion is that the code is working.
Just to re-iterate, the atoms to be pulled to the bilayer center were present at the lipid-water interface in both leaflets. Here is the (apparently successful) pull code: ; COM PULLING pull = umbrella pull_geometry = position pull_dim = N N Y pull_start = no pull_nstxout = 5000 pull_nstfout = 5000 pull_ngroups = 2 pull_group0 = POPC ; global atom number of a POPC atom in the middle of the xy box pull_pbcatom0 = 3100 pull_group1 = C35TOP ; global atom number of a C35 atom in the middle of the xy box pull_pbcatom1 = 5467 pull_init1 = 0 0 0.0 pull_rate1 = 0 pull_k1 = 500 pull_vec1 = 0 0 0 pull_group2 = C35BOT ; global atom number of a C35 atom in the middle of the xy box pull_pbcatom2 = 7755 pull_init2 = 0 0 0.0 pull_rate2 = 0 pull_k2 = 500 pull_vec2 = 0 0 0 ########################## Here is the output from grompp Pull group natoms pbc atom distance at start reference at t=0 0 4992 3100 1 16 5467 0.000 0.000 3.184 0.000 0.000 0.000 2 16 7755 -0.000 0.000 -3.160 0.000 0.000 0.000 ########################### The average COM pull En. and potential from the log file after 100 ps. 2.96054e+03 ########################### And here is the output in the log files, detailing the pulling parameters pull = umbrella pull_geometry = position pull_dim (3): pull_dim[0]=0 pull_dim[1]=0 pull_dim[2]=1 pull_r1 = 1 pull_r0 = 1.5 pull_constr_tol = 1e-06 pull_nstxout = 5000 pull_nstfout = 5000 pull_ngrp = 2 pull_group 0: atom (4992): atom[0,...,4991] = {0,...,4991} weight: not available pbcatom = 3099 vec (3): vec[0]= 0.00000e+00 vec[1]= 0.00000e+00 vec[2]= 0.00000e+00 init (3): init[0]= 0.00000e+00 init[1]= 0.00000e+00 init[2]= 0.00000e+00 rate = 0 k = 0 kB = 0 pull_group 1: atom (16): atom[0]=5026 atom[1]=5114 atom[2]=5202 atom[3]=5290 atom[4]=5378 atom[5]=5466 atom[6]=5554 atom[7]=5642 atom[8]=5730 atom[9]=5818 atom[10]=5906 atom[11]=5994 atom[12]=6082 atom[13]=6169 atom[14]=6258 atom[15]=6346 weight: not available pbcatom = 5466 vec (3): vec[0]= 0.00000e+00 vec[1]= 0.00000e+00 vec[2]= 0.00000e+00 init (3): init[0]= 0.00000e+00 init[1]= 0.00000e+00 init[2]= 0.00000e+00 rate = 0 k = 500 kB = 500 pull_group 2: atom (16): atom[0]=6434 atom[1]=6522 atom[2]=6610 atom[3]=6698 atom[4]=6786 atom[5]=6874 atom[6]=6962 atom[7]=7050 atom[8]=7138 atom[9]=7226 atom[10]=7314 atom[11]=7402 atom[12]=7490 atom[13]=7578 atom[14]=7666 atom[15]=7754 weight: not available pbcatom = 7754 vec (3): vec[0]= 0.00000e+00 vec[1]= 0.00000e+00 vec[2]= 0.00000e+00 init (3): init[0]= 0.00000e+00 init[1]= 0.00000e+00 init[2]= 0.00000e+00 rate = 0 k = 500 kB = 500 ######################### I would want to make sure I am doing this right Best Wishes Maria On Wed, Sep 8, 2010 at 4:09 PM, <chris.ne...@utoronto.ca> wrote: > Maria, > > I highly doubt that you are correct. If you want more help, please copy and > paste your commands and some relevant output back to the list. > > Chris. > > --original message -- > > Message: 4 > Date: Wed, 8 Sep 2010 15:07:15 +0200 > From: maria goranovic <mariagorano...@gmail.com> > Subject: Re: [gmx-users] restraining atoms to the plane at the bilayer > center: pull code ? > To: Discussion list for GROMACS users <gmx-users@gromacs.org> > Message-ID: > <aanlkti=rqv+ozsff0u5pmxsuaq-jpjaghnnvpdrjm...@mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > Dear Chris, > > Thank you for the code. I did check out pull_pbcatomN, and thank you for > the > heads up that it is a global index. I used a central atom in the box for > each of the two groups as pull_pbcatom. > > After using the code below with a rate constant of ~ 5500 for about 100 ps, > I realized that the atoms in different leaflets were being pulled in the > same direction. I wanted certain atoms of either leaflet to be restrained > to the bilayer center. However, the code pulls all the atoms in the same > direction, such that atoms from one leaflet do approach the bilayer center, > but those from the other leaflet tend to get further away from the bilayer > center. Thus, one actually needs two groups to be pulled, one each for the > atoms in each leaflet. Once this was fixed, the code worked just fine. > > Again, thank you. > > Maria > > > > On Tue, Sep 7, 2010 at 3:40 PM, <chris.ne...@utoronto.ca> wrote: > > [Hide Quoted Text] > Maria, try this. There actually is a lot of this on the mailing list, so I > suggest checking it a little deeper for your next querry, or at least > outlining how you looked and what you found so that it is clear you have > tried. > > Also, read about pull_pbcatomN and think carefully about how you want to > set that up. It is a *global* index. > > ; COM PULLING > pull = umbrella > pull_geometry = position > pull_dim = N N Y > pull_start = no > pull_nstxout = 500 > pull_nstfout = 500 > pull_ngroups = 1 > pull_group0 = <<<bilayer-group>>> > pull_pbcatom0 = 0 > pull_group1 = <<<pulled-group>>> > pull_pbcatom1 = 0 > pull_init1 = 0 0 0.0 > pull_rate1 = 0 > pull_k1 = <<<force-constant>>> > pull_vec1 = 0 0 0 > > -- original message -- > > > The manual does discuss restraining to a plane, but this must be the plane > in which the atom is already present. > > [ position_restraints ] > ; ai funct fc > ... > 3 1 1000 0 0 > > How about restraining the atom to some other plane? For example, how about > restraining a phosphate group initially at the lipid-water interface to the > bilayer center (for whatever fancy reasons) ? Won't this require pulling > the > atom to that plane first? > > If it can be achieved using the position restraints alone, it is not clear > to me how to do this? > > -Maria > > > > On Tue, Sep 7, 2010 at 12:53 PM, Justin A. Lemkul <jalemkul at vt.edu> > wrote: > > maria goranovic wrote: > > I want to restrain certain atoms of my simulation to the plane > perpendicular to the bilayer normal, and at the bilayer center. Can > someone > please provide a quick guide on how to do this? I read the pull-code > options, but restraining to a plane did not seem possible? > > > You can restrain atoms to a plane using position restraints. See the > manual for an example. > -Justin > -- > Maria G. > Technical University of Denmark > Copenhagen > > -- > gmx-users mailing list gmx-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 thewww interface > or send it to gmx-users-requ...@gromacs.org. > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > -- Maria G. Technical University of Denmark Copenhagen
-- gmx-users mailing list gmx-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