<html>
<head>
<style>
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Verdana
}
</style>
</head>
<body class='hmmessage'>
Hi,<br><br>Just to be sure, the pull code is producing the correct results now?<br>The only problems should be in several analysis tools.<br><br>The analysis tools have been written only with analysis of molecules<br>or parts of molecules in mind. When you want to analyze groups<br>that cover two molecules or more there will be pbc problems.<br>These are not simple to fix. If the distances between all the atoms<br>in the group are less than half a box length there is a unique COM.<br>The procedure should be similar to what trjconv -pbc cluster does.<br><br>I can try to implement this for g_traj and g_dist.<br>Are there other tools that cause problems?<br><br>Berk<br><br>> Date: Mon, 26 Jan 2009 17:46:26 -0500<br>> From: fiedler@umich.edu<br>> To: gmx-users@gromacs.org<br>> Subject: Re: [gmx-users] Transition difficulties: version 3.3.3 to 4.0.3 regarding pull_geometry=distance<br>> <br>> Hi Chris,<br>> <br>> Your suggestions were helpful. Due to the nature of this problem, there <br>> were limitations in the application of the sequence of steps that you <br>> recommended. Specifically, the difference in the equilibrated structure <br>> from version 3.3.x to the calculated reference center of mass from <br>> version 4.0.x causes an unreasonable "jump" with constraint force <br>> calculation. Resultant bad contacts prohibit generation of a pull <br>> output file.<br>> <br>> More generally though, I believe I understand your broader point: the <br>> sensitive nature of the periodic box with respect to the pull code and <br>> possible discrepancies in the pbc treatment by various utilities, may <br>> cause confusion with the diagnostic process. The test bilayer system I <br>> created however, centered distantly from the x and y box edges, may <br>> effectively remove this set of concerns. For example, this <br>> configuration is relatively unchanged upon "recentering", however it is <br>> still susceptible to the problem that is the subject of this discussion.<br>> <br>> If I can reproduce this phenomenon with a publicly available topology <br>> file, I can provide the necessary input files for inspection.<br>> <br>> Thank you again,<br>> <br>> Steve Fiedler<br>> <br>> <br>> <br>> chris.neale@utoronto.ca wrote:<br>> > Hi Steve,<br>> ><br>> > what I intended to suggest was actually something different (and much <br>> > easier).<br>> ><br>> > The idea is not that you need some special system to be able to <br>> > utilize the pull code, but that the pull code is correct whereas the <br>> > g_dist and g_traj programs are not as good at treating pbc in the way <br>> > that one desires.<br>> ><br>> > I suggest the following.<br>> ><br>> > 1. Take your original system and run the pull code for a very short <br>> > simulation. Use the last line of the output to calculate the relevant <br>> > displacement<br>> ><br>> > 2. Now use trjconv -b -e to get the last frame of the .xtc that <br>> > resulted from that short MD run as a .gro file, call it final.gro. I <br>> > suspect that your groups are not entirely in the same simulation box <br>> > in final.gro.<br>> ><br>> > 3. Now make a new .ndx file from that .gro and give it a single <br>> > residue that is near your binding pocket, call it R_1<br>> ><br>> > 4. Now apply trjconv -center -pbc mol -ur compact while selecting R_1 <br>> > for centering, call the new .gro file final_center.gro<br>> ><br>> > 5. Visualize final_center.gro and ensure that all of your relevant <br>> > atoms are in the same image in the way that puts the minimum distance <br>> > between them along a path that is entirely contained within the unit <br>> > cell. If not, go back to step 3 and try making a group R_2, etc. <br>> > until this process works. NOTE: you might think that giving trjconv <br>> > -center the relevant groups that you use for pulling will be a good <br>> > idea here, but it is not. The problem there is that the atoms may be <br>> > "centered" by placing half on the left boundary and half on the right <br>> > boundary. I find using one logically selected residue or atom is the <br>> > best method here.<br>> ><br>> > 6. Assuming that you got what you wanted in step 5, now run g_traj and <br>> > g_dist on final_center.gro. In my case, I found that g_traj and g_dist <br>> > give the same answer as the pull code output when I am using <br>> > final_center.gro, but not always when I am using final.gro.<br>> ><br>> > *** I always laugh when these problems arise because, in an important <br>> > sense, the protein *did* jump out of the simulation box... at least as <br>> > far as g_traj and g_dist are concerned. This, we must hope, is <br>> > correctly treated in the pull code even though it is incorrectly (or <br>> > at least unintuitively) treated by g_traj and g_dist.<br>> ><br>> > Chris.<br>> ><br>> > -- original message --<br>> ><br>> > Hi,<br>> ><br>> > Thank you Berk and Chris for the suggestions.<br>> ><br>> > To address the possibility that this issue is related to periodic<br>> > boundaries, I used two approaches:<br>> > 1. The pull group of interest (permeant) was centered in the x-y plane<br>> > of the box using Chris' approach. I then used the genconf utility to<br>> > replicate my lipid box to a 9x9 grid in the x-y plane and removed all<br>> > but the center box. This generated the coordinates for a bilayer system<br>> > with all lipid molecules inside a box and intact. The discrepancy<br>> > between the grompp (version 4.0.3) output and distances as calculated by<br>> > g_traj (version 4.0.3) persist, 2.667 vs. 0.3996 nm.<br>> > 2. I constructed a three atom system containing 2 reference atoms of<br>> > type A, and a "pull" atom of type B. Proper output from grompp was<br>> > observed for all coordinates of both the reference and pulled atoms,<br>> > include coordinates for atoms moved outside the box in the x-y plane.<br>> > The coordinate, topology, and run control parameter file are given below.<br>> ><br>> > If there are additional suggestions, I would be greatly appreciative.<br>> ><br>> > Thank you,<br>> ><br>> > Steve Fiedler<br>> ><br>> > -----------------<br>> > conf.gro<br>> > Three atoms<br>> > 3<br>> > 1AAA A 1 1.500 1.500 1.000<br>> > 2AAA A 2 0.500 1.500 1.000<br>> > 3BBB B 3 -1.500 1.500 1.700<br>> > 3.00000 3.00000 3.00000<br>> > -----------------<br>> > index.ndx<br>> > [ System ]<br>> > 1 2 3<br>> > [ Ref ]<br>> > 1 2<br>> > [ Pulled ]<br>> > 3<br>> > -----------------<br>> > grompp.mdp<br>> > title = ThreeAtoms<br>> > integrator = md<br>> > dt = 0.001<br>> > nsteps = 1<br>> > ns_type = grid<br>> > pbc = xyz<br>> > coulombtype = shift<br>> > rlist = 1.4<br>> > rcoulomb = 1.4<br>> > rvdw = 1.4<br>> > tcoupl = no<br>> > pcoupl = no<br>> > constraint_algorithm = shake<br>> > shake_tol = 1e-4<br>> > gen-vel = no<br>> > gen-temp = 0<br>> ><br>> > nstxout = 1<br>> > nstvout = 0<br>> > nstfout = 0<br>> ><br>> > pull = umbrella<br>> > pull_geometry = distance<br>> > pull_dim = N N Y<br>> > pull_start = no<br>> > pull_init1 = 0.7<br>> > pull_group0 = Ref<br>> > pull_group1 = Pulled<br>> > pull_k1 = 10000<br>> > -----------------<br>> > topology.top<br>> > ; topology for two partially charged atoms<br>> ><br>> > [ defaults ]<br>> > ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ<br>> > 1 3 yes 0.125 0.5<br>> ><br>> > [ atomtypes ]<br>> > ;name mass charge ptype sig eps<br>> > A 1000.0000 0.000 A 0.50000 9.90000<br>> > B 9.0000 0.000 A 0.30000 9.00000<br>> ><br>> > [ nonbond_params ]<br>> > ; i j func sig eps<br>> ><br>> > [ moleculetype ]<br>> > AAAA 1<br>> ><br>> > [ atoms ]<br>> > ; nr type resnr residue atom cgnr charge mass<br>> > 1 A 1 AAA A 1 0.000 1000.0000<br>> ><br>> > [ moleculetype ]<br>> > BBBB 1<br>> ><br>> > [ atoms ]<br>> > ; nr type resnr residue atom cgnr charge mass<br>> > 1 B 1 BBB B 1 0.000 9.00<br>> ><br>> > [ system ]<br>> > ; name<br>> > Three atoms<br>> ><br>> > [ molecules ]<br>> > ; name number<br>> > AAAA 2<br>> > BBBB 1<br>> ><br>> ><br>> ><br>> ><br>> ><br>> > Chris Neale wrote:<br>> >> I just checked similar simulations of mine and Berk's suggestion <br>> >> accounts for similar discrepancies that I notice on a quick <br>> >> evaluation where g_traj and g_dist fail to give me the same distance <br>> >> as I obtain from the pull pos.xvg file. As Berk suggests, once I <br>> >> first trjconv -center -pbc mol -ur compact (giving an appropriate <br>> >> residue for centering that puts all relevant pulled atoms in the same <br>> >> box) then g_traj and g_dist both give me the exact same answer as I <br>> >> calculate based on pull pos.xvg. Chris -- original message -- Hi, <br>> >> There could be a problem with periodic boundary conditions. Do you <br>> >> have multiple molecules in a pull group, or broken molecules? In that <br>> >> case the COM position of 3.3.3 and g_traj are both incorrect. The <br>> >> pull code in 4.0 grompp and mdrun are (as far as I know) always <br>> >> correct. Berk<br>> ><br>> > _______________________________________________<br>> > gmx-users mailing list gmx-users@gromacs.org<br>> > http://www.gromacs.org/mailman/listinfo/gmx-users<br>> > Please search the archive at http://www.gromacs.org/search before <br>> > posting!<br>> > Please don't post (un)subscribe requests to the list. Use thewww <br>> > interface or send it to gmx-users-request@gromacs.org.<br>> > Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br>> ><br>> ><br>> <br>> <br>> -- <br>> Steve Fiedler, Ph.D.<br>> Research Fellow<br>> Department of Mechanical Engineering<br>> The University of Michigan<br>> 2024 G.G. Brown<br>> 2350 Hayward St.<br>> Ann Arbor, MI 48109-2125<br>> <br>> <br>> _______________________________________________<br>> gmx-users mailing list gmx-users@gromacs.org<br>> http://www.gromacs.org/mailman/listinfo/gmx-users<br>> Please search the archive at http://www.gromacs.org/search before posting!<br>> Please don't post (un)subscribe requests to the list. Use the <br>> www interface or send it to gmx-users-request@gromacs.org.<br>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br><br /><hr />What can you do with the new Windows Live? <a href='http://www.microsoft.com/windows/windowslive/default.aspx' target='_new'>Find out</a></body>
</html>