[gmx-users] Re: a flat plane boundary condition
David
spoel at xray.bmc.uu.se
Wed May 25 23:16:49 CEST 2005
On Wed, 2005-05-25 at 16:59 -0400, Beckett W Sterner wrote:
> Hi,
> Specifically, the biology is such that the C terminus is exiting a pore
> in the cell's external membrane, and in vivo the protein folds from the C
> terminus as it exits. The process of the protein exiting the pore is
> about 30-40 residues per second, however, for a protein of >200 residues.
> Since the process is so slow, I want to treat the part of the chain still
> inside the pore as essentially fixed in place (and hence not worth
> simulating). So the idea for the simulation is that if there are X
> residues that have exited the pore, I want to let that part fold for a
> short period of time (<< 1/40 seconds), and then plug that result into
> another run with X+Y residues, so that Y residues have effectively
> "exited" the pore between the X and X+Y runs. I would repeat this until X
> equaled a substantial portion of the protein.
Since you're not modeling the pore it self, it seems that is not very
important what is there. You could e.g. build a (frozen or position
restrained) hydrophobic wall to which you attach the N-terminal end of
the Y peptide.
>
> Thanks,
> Beckett
>
>
> On Wed, 25 May 2005, David wrote:
>
> > On Wed, 2005-05-25 at 15:56 -0400, Beckett W Sterner wrote:
> > > Hi,
> > > Thats looks great for fixing one end of the protein, but I also want to
> > > ensure that the rest of it is restricted to only half of space (the
> > > impenetrable plane boundary condition). Or maybe I'm missing your
> > > meaning. The constraint is important because as the protein exits the
> > > pore it should begin to fold but it shouldn't clump around the residue
> > > with fixed position.
> > >
> > It's not clear to me what you want, or at least how you model the
> > protein coming out of the pore. Do you want to allow an ever increasing
> > (N-terminal) end to start folding in solution while you restrain the C-
> > terminal end?
> >
> > > Thanks,
> > > Beckett
> > >
> > > > You can freeze atoms in one, two or three dimensions. Check mdp options
> > > > on www.gromacs.org
> > >
> > >
> > >
> > >
> > > > On Wed, 2005-05-25 at 12:30 -0400, Beckett W Sterner wrote:
> > > > > Hi,
> > > > > I'm interested in doing a simple simulation of how a protein folds as it
> > > > > slowly emerges from a pore through a cell membrane. I unfortunately do
> > > > > not have a crystal structure for the pore, and would like to keep
> > > > > computational costs to a minimum, so I'm interested in running an MD
> > > > > simulation on the protein where one end is fixed to a flat plane with a 2d
> > > > > periodic boundary condition. Is there a way to implement this in GROMACS,
> > > > > perhaps using a force field?
> > >
> > >
> > >
> > > > You can freeze atoms in one, two or three dimensions. Check mdp options
> > > > on www.gromacs.org
> > > > >
> > > > > Thanks,
> > > > > Beckett
> > > > > _______________________________________________
> > > > > gmx-users mailing list
> > > > > gmx-users at gromacs.org
> > > > > http://www.gromacs.org/mailman/listinfo/gmx-users
> > > > > Please don't post (un)subscribe requests to the list. Use the
> > > > > www interface or send it to gmx-users-request at gromacs.org.
> > > > --
> > > > David.
> > > > ________________________________________________________________________
> > > > David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
> > > > Dept. of Cell and Molecular Biology, Uppsala University.
> > > > Husargatan 3, Box 596, 75124 Uppsala, Sweden
> > > > phone: 46 18 471 4205 fax: 46 18 511 755
> > > > spoel at xray.bmc.uu.se spoel at gromacs.org http://xray.bmc.uu.se/~spoel
> > > > ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> > > >
> > > >
> > > >
> > > > ------------------------------
> > > >
> > > > Message: 3
> > > > Date: Wed, 25 May 2005 13:55:46 -0400
> > > > From: Anthony Cruz <acb15885 at uprm.edu>
> > > > Subject: Re: [gmx-users] g_rms Fatal error: Error: Too many iterations
> > > > in routine JACOBI
> > > > To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> > > > Message-ID: <200505251355.46962.acb15885 at uprm.edu>
> > > > Content-Type: text/plain; charset="iso-8859-1"
> > > >
> > > > I do the gmxcheck and this is the result:
> > > >
> > > > Checking file 1scn_bx_w_fullMD_mv.trr
> > > > trn version: GMX_trn_file (single precision)
> > > > Reading frame 0 time 0.000
> > > > # Atoms 110719
> > > > Last frame 10500 time 10500.001
> > > >
> > > >
> > > > Item #frames Timestep (ps)
> > > > Step 10501 1
> > > > Time 10501 1
> > > > Lambda 10501 1
> > > > Coords 10501 1
> > > > Velocities 10501 1
> > > > Forces 0
> > > > Box 10501 1
> > > > Checking coordinate file 1scn_bx_w_CL_TPR4fullMD.tpr
> > > > Reading file 1scn_bx_w_CL_TPR4fullMD.tpr, VERSION 3.2.1 (single precision)
> > > > Reading file 1scn_bx_w_CL_TPR4fullMD.tpr, VERSION 3.2.1 (single precision)
> > > > 110719 atoms in file
> > > > coordinates found
> > > > box found
> > > > velocities found
> > > >
> > > > Kinetic energy: 411492 (kJ/mol)
> > > > Assuming the number of degrees of freedom to be Natoms * 3 or Natoms * 2,
> > > > the velocities correspond to a temperature of the system
> > > > of 297.996 K or 446.995 K respectively.
> > > >
> > > > Checking for atoms closer than 0.8 and not between 0.4 and 0.7,
> > > > relative to sum of Van der Waals distance:
> > > > WARNING: masses will be determined based on residue and atom names,
> > > > this can deviate from the real mass of the atom type
> > > > In case you use free energy of solvation predictions:
> > > >
> > > > ++++++++ PLEASE CITE THE FOLLOWING REFERENCE ++++++++
> > > > D. Eisenberg and A. D. McLachlan
> > > > Solvation energy in protein folding and binding
> > > > Nature 319 (1986) pp. 199-203
> > > > -------- -------- --- Thank You --- -------- --------
> > > >
> > > > Opening library file /usr/local/gromacs/share/top/aminoacids.dat
> > > > Opening library file /usr/local/gromacs/share/top/atommass.dat
> > > > Opening library file /usr/local/gromacs/share/top/vdwradii.dat
> > > > Opening library file /usr/local/gromacs/share/top/dgsolv.dat
> > > > #Entries in atommass.dat: 82 vdwradii.dat: 26 dgsolv.dat: 7
> > > > atom# name residue r_vdw atom# name residue r_vdw distance
> > > > 40 CB PRO 5 0.15 42 CD PRO 5 0.15 0.2332
> > > > 50 CD1 TYR 6 0.15 58 CZ TYR 6 0.15 0.2369
> > > > 52 CD2 TYR 6 0.15 58 CZ TYR 6 0.15 0.2358
> > > > 78 CA PRO 9 0.15 80 CG PRO 9 0.15 0.2383
> > > > 87 CB LEU 10 0.15 91 C LEU 10 0.15 0.2399
> > > > 149 C VAL 16 0.15 153 CA GLN 17 0.15 0.2398
> > > > 190 CG PHE 21 0.15 195 CE1 PHE 21 0.15 0.2395
> > > > 191 CD1 PHE 21 0.15 193 CD2 PHE 21 0.15 0.2393
> > > > 193 CD2 PHE 21 0.15 199 CZ PHE 21 0.15 0.2369
> > > > 244 C VAL 26 0.15 248 CA LYSH 27 0.15 0.2371
> > > > 277 CG1 VAL 30 0.15 278 CG2 VAL 30 0.15 0.2398
> > > > 302 CB THR 33 0.15 306 C THR 33 0.15 0.2383
> > > > 352 CG HISB 39 0.15 355 CE1 HISB 39 0.15 0.2174
> > > > 354 CD2 HISB 39 0.15 355 CE1 HISB 39 0.15 0.2169
> > > > 361 CA PRO 40 0.15 364 CD PRO 40 0.15 0.2395
> > > > 362 CB PRO 40 0.15 364 CD PRO 40 0.15 0.2387
> > > > 383 C LEU 42 0.15 387 CA ASN 43 0.15 0.2391
> > > > 440 CG PHE 50 0.15 447 CE2 PHE 50 0.15 0.2371
> > > > 441 CD1 PHE 50 0.15 449 CZ PHE 50 0.15 0.2331
> > > > 443 CD2 PHE 50 0.15 449 CZ PHE 50 0.15 0.2371
> > > > 459 C VAL 51 0.15 463 CA ALA 52 0.15 0.2376
> > > > 492 CG TYR 56 0.15 497 CE1 TYR 56 0.15 0.2334
> > > > 492 CG TYR 56 0.15 499 CE2 TYR 56 0.15 0.2398
> > > > 497 CE1 TYR 56 0.15 499 CE2 TYR 56 0.15 0.2375
> > > > 558 CA HISA 63 0.15 560 CG HISA 63 0.15 0.2398
> > > > 560 CG HISA 63 0.15 564 CE1 HISA 63 0.15 0.212
> > > > 563 CD2 HISA 63 0.15 564 CE1 HISA 63 0.15 0.2126
> > > > 586 CG HISB 66 0.15 589 CE1 HISB 66 0.15 0.2163
> > > > 588 CD2 HISB 66 0.15 589 CE1 HISB 66 0.15 0.2109
> > > > 692 C GLY 79 0.15 696 CA VAL 80 0.15 0.2391
> > > > 698 CG1 VAL 80 0.15 699 CG2 VAL 80 0.15 0.2384
> > > > 732 CB PRO 85 0.15 734 CD PRO 85 0.15 0.2341
> > > > 751 C VAL 87 0.15 755 CA SER 88 0.15 0.2397
> > > > 775 CD1 TYR 90 0.15 783 CZ TYR 90 0.15 0.2391
> > > > 777 CD2 TYR 90 0.15 783 CZ TYR 90 0.15 0.2398
> > > > 817 CA VAL 94 0.15 819 CG1 VAL 94 0.15 0.2397
> > > > 890 CD1 TYR 103 0.15 892 CD2 TYR 103 0.15 0.2384
> > > > 894 CE1 TYR 103 0.15 896 CE2 TYR 103 0.15 0.2346
> > > > 929 CG1 VAL 107 0.15 930 CG2 VAL 107 0.15 0.2398
> > > > 969 CG TRP 112 0.15 975 CE2 TRP 112 0.15 0.2231
> > > > 970 CD1 TRP 112 0.15 972 CD2 TRP 112 0.15 0.2171
> > > > 970 CD1 TRP 112 0.15 975 CE2 TRP 112 0.15 0.2157
> > > > 975 CE2 TRP 112 0.15 976 CE3 TRP 112 0.15 0.2367
> > > > 978 CZ2 TRP 112 0.15 980 CZ3 TRP 112 0.15 0.2378
> > > > 1079 C MET 123 0.15 1083 CA SER 124 0.15 0.2396
> > > > 1094 CD1 LEU 125 0.15 1095 CD2 LEU 125 0.15 0.239
> > > > 1216 C ASN 140 0.15 1220 CA ALA 141 0.15 0.2389
> > > > 1228 CG TYR 142 0.15 1233 CE1 TYR 142 0.15 0.2379
> > > > 1229 CD1 TYR 142 0.15 1231 CD2 TYR 142 0.15 0.2377
> > > > 1233 CE1 TYR 142 0.15 1235 CE2 TYR 142 0.15 0.2356
> > > > 1251 CB ARG 144 0.15 1263 C ARG 144 0.15 0.2347
> > > > 1347 C GLY 156 0.15 1351 CA ASN 157 0.15 0.2398
> > > > 1428 CG TYR 166 0.15 1433 CE1 TYR 166 0.15 0.238
> > > > 1428 CG TYR 166 0.15 1435 CE2 TYR 166 0.15 0.2363
> > > > 1443 CA PRO 167 0.15 1446 CD PRO 167 0.15 0.2391
> > > > 1473 CD1 TYR 170 0.15 1475 CD2 TYR 170 0.15 0.237
> > > > 1473 CD1 TYR 170 0.15 1481 CZ TYR 170 0.15 0.2383
> > > > 1477 CE1 TYR 170 0.15 1479 CE2 TYR 170 0.15 0.2392
> > > > 1530 CG1 VAL 176 0.15 1531 CG2 VAL 176 0.15 0.2373
> > > > 1635 CG PHE 188 0.15 1640 CE1 PHE 188 0.15 0.2336
> > > > 1635 CG PHE 188 0.15 1642 CE2 PHE 188 0.15 0.2393
> > > > 1640 CE1 PHE 188 0.15 1642 CE2 PHE 188 0.15 0.2399
> > > > 1686 CB GLU 194 0.15 1691 C GLU 194 0.15 0.2396
> > > > 1698 CD1 LEU 195 0.15 1699 CD2 LEU 195 0.15 0.2293
> > > > 1736 CA PRO 200 0.15 1738 CG PRO 200 0.15 0.2346
> > > > 1762 CG1 VAL 204 0.15 1763 CG2 VAL 204 0.15 0.2344
> > > > 1771 CD1 TYR 205 0.15 1773 CD2 TYR 205 0.15 0.2366
> > > > 1775 CE1 TYR 205 0.15 1777 CE2 TYR 205 0.15 0.2391
> > > > 1806 CD1 TYR 208 0.15 1808 CD2 TYR 208 0.15 0.2383
> > > > 1810 CE1 TYR 208 0.15 1812 CE2 TYR 208 0.15 0.2356
> > > > 1859 CG TYR 213 0.15 1864 CE1 TYR 213 0.15 0.2375
> > > > 1860 CD1 TYR 213 0.15 1862 CD2 TYR 213 0.15 0.24
> > > > 1862 CD2 TYR 213 0.15 1868 CZ TYR 213 0.15 0.2396
> > > > 1864 CE1 TYR 213 0.15 1866 CE2 TYR 213 0.15 0.2366
> > > > 1954 CA PRO 224 0.15 1957 CD PRO 224 0.15 0.2351
> > > > 1964 CG HISA 225 0.15 1968 CE1 HISA 225 0.15 0.2137
> > > > 1967 CD2 HISA 225 0.15 1968 CE1 HISA 225 0.15 0.2127
> > > > 2012 CB LEU 232 0.15 2014 CD1 LEU 232 0.15 0.2396
> > > > 2032 CD1 LEU 234 0.15 2033 CD2 LEU 234 0.15 0.2395
> > > > 2063 CD2 HISB 237 0.15 2064 CE1 HISB 237 0.15 0.2173
> > > > 2070 CA PRO 238 0.15 2072 CG PRO 238 0.15 0.2387
> > > > 2236 CG TYR 255 0.15 2243 CE2 TYR 255 0.15 0.2388
> > > > 2241 CE1 TYR 255 0.15 2243 CE2 TYR 255 0.15 0.237
> > > > 2284 CG PHE 260 0.15 2289 CE1 PHE 260 0.15 0.2363
> > > > 2285 CD1 PHE 260 0.15 2287 CD2 PHE 260 0.15 0.2389
> > > > 2287 CD2 PHE 260 0.15 2293 CZ PHE 260 0.15 0.2339
> > > > 2301 CG TYR 261 0.15 2308 CE2 TYR 261 0.15 0.2389
> > > > 2302 CD1 TYR 261 0.15 2310 CZ TYR 261 0.15 0.2358
> > > > 2304 CD2 TYR 261 0.15 2310 CZ TYR 261 0.15 0.2395
> > > > 2319 CG TYR 262 0.15 2326 CE2 TYR 262 0.15 0.2376
> > > > 2324 CE1 TYR 262 0.15 2326 CE2 TYR 262 0.15 0.238
> > > >
> > > > Atoms outside box ( 10.4513 10.4513 10.4513 ):
> > > > (These may occur often and are normally not a problem)
> > > > atom# name residue r_vdw coordinate
> > > > 2488 OW SOL 293 0.105 -0.034 2.79 4.15
> > > > 2491 OW SOL 294 0.105 -0.031 5.54 9.37
> > > > 2494 OW SOL 295 0.105 -0.024 3.92 10.1
> > > > 2497 OW SOL 296 0.105 -0.02 5.24 10.3
> > > > 2500 OW SOL 297 0.105 -0.018 4.16 8.09
> > > > 2501 HW1 SOL 297 0.04 -0.027 4.18 8.19
> > > > 2503 OW SOL 298 0.105 -0.018 6.61 5.17
> > > > 2506 OW SOL 299 0.105 -0.016 1.37 7.5
> > > > 2509 OW SOL 300 0.105 -0.015 2.23 9.93
> > > > 2511 HW2 SOL 300 0.04 -0.048 2.29 10
> > > > (maybe more)-bash-2.05b$ g_rms -f 1scn_bx_w_fullMD_mv.trr -s
> > > > 1scn_bx_w_CL_TPR4fullMD.tpr -o
> > > > :-) G R O M A C S (-:
> > > >
> > > > Grunge ROck MAChoS
> > > >
> > > > :-) VERSION 3.2.1 (-:
> > > >
> > > >
> > > > Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
> > > > Copyright (c) 1991-2000, University of Groningen, The Netherlands.
> > > > Copyright (c) 2001-2004, The GROMACS development team,
> > > > check out http://www.gromacs.org for more information.
> > > >
> > > > This program is free software; you can redistribute it and/or
> > > > modify it under the terms of the GNU General Public License
> > > > as published by the Free Software Foundation; either version 2
> > > > of the License, or (at your option) any later version.
> > > >
> > > > :-) g_rms (-:
> > > >
> > > > Option Filename Type Description
> > > > ------------------------------------------------------------
> > > > -s 1scn_bx_w_CL_TPR4fullMD.tpr Input Structure+mass(db): tpr tpb
> > > > tpa gro g96 pdb xml
> > > > -f 1scn_bx_w_fullMD_mv.trr Input Generic trajectory: xtc trr trj
> > > > gro g96 pdb
> > > > -f2 traj.xtc Input, Opt. Generic trajectory: xtc trr trj gro g96 pdb
> > > > -n index.ndx Input, Opt. Index file
> > > > -o rmsd.xvg Output xvgr/xmgr file
> > > > -mir rmsdmir.xvg Output, Opt. xvgr/xmgr file
> > > > -a avgrp.xvg Output, Opt. xvgr/xmgr file
> > > > -dist rmsd-dist.xvg Output, Opt. xvgr/xmgr file
> > > > -m rmsd.xpm Output, Opt. X PixMap compatible matrix file
> > > > -bin rmsd.dat Output, Opt. Generic data file
> > > > -bm bond.xpm Output, Opt. X PixMap compatible matrix file
> > > >
> > > > Option Type Value Description
> > > > ------------------------------------------------------
> > > > -[no]h bool no Print help info and quit
> > > > -[no]X bool no Use dialog box GUI to edit command line options
> > > > -nice int 19 Set the nicelevel
> > > > -b time -1 First frame (ps) to read from trajectory
> > > > -e time -1 Last frame (ps) to read from trajectory
> > > > -dt time -1 Only use frame when t MOD dt = first time (ps)
> > > > -tu enum ps Time unit: ps, fs, ns, us, ms, s, m or h
> > > > -[no]w bool no View output xvg, xpm, eps and pdb files
> > > > -what enum rmsd Structural difference measure: rmsd, rho or rhosc
> > > > -[no]pbc bool yes PBC check
> > > > -fit enum rot+trans Fit to reference structure: rot+trans,
> > > > translation or none
> > > > -prev int 0 Compare with previous frame
> > > > -[no]split bool no Split graph where time is zero
> > > > -skip int 1 Only write every nr-th frame to matrix
> > > > -skip2 int 1 Only write every nr-th frame to matrix
> > > > -max real -1 Maximum level in comparison matrix
> > > > -min real -1 Minimum level in comparison matrix
> > > > -bmax real -1 Maximum level in bond angle matrix
> > > > -bmin real -1 Minimum level in bond angle matrix
> > > > -nlevels int 80 Number of levels in the matrices
> > > >
> > > > Reading file 1scn_bx_w_CL_TPR4fullMD.tpr, VERSION 3.2.1 (single precision)
> > > > Reading file 1scn_bx_w_CL_TPR4fullMD.tpr, VERSION 3.2.1 (single precision)
> > > > Select group for least squares fit
> > > > Opening library file /usr/local/gromacs/share/top/aminoacids.dat
> > > > Group 0 ( System) has 110719 elements
> > > > Group 1 ( Protein) has 2433 elements
> > > > Group 2 ( Protein-H) has 1920 elements
> > > > Group 3 ( C-alpha) has 274 elements
> > > > Group 4 ( Backbone) has 822 elements
> > > > Group 5 ( MainChain) has 1097 elements
> > > > Group 6 (MainChain+Cb) has 1336 elements
> > > > Group 7 ( MainChain+H) has 1364 elements
> > > > Group 8 ( SideChain) has 1069 elements
> > > > Group 9 ( SideChain-H) has 823 elements
> > > > Group 10 ( Prot-Masses) has 2433 elements
> > > > Group 11 ( Non-Protein) has 108286 elements
> > > > Group 12 ( CL) has 4 elements
> > > > Group 13 ( HOH) has 423 elements
> > > > Group 14 ( SOL) has 107856 elements
> > > > Group 15 ( CA) has 2 elements
> > > > Group 16 ( NA) has 1 elements
> > > > Group 17 ( Other) has 108286 elements
> > > > Select a group: 1
> > > > Selected 1: 'Protein'
> > > > How many groups do you want to compare ? 1
> > > > OK. I will compare 1 group
> > > > Select group for RMSD calculation
> > > > Opening library file /usr/local/gromacs/share/top/aminoacids.dat
> > > > Group 0 ( System) has 110719 elements
> > > > Group 1 ( Protein) has 2433 elements
> > > > Group 2 ( Protein-H) has 1920 elements
> > > > Group 3 ( C-alpha) has 274 elements
> > > > Group 4 ( Backbone) has 822 elements
> > > > Group 5 ( MainChain) has 1097 elements
> > > > Group 6 (MainChain+Cb) has 1336 elements
> > > > Group 7 ( MainChain+H) has 1364 elements
> > > > Group 8 ( SideChain) has 1069 elements
> > > > Group 9 ( SideChain-H) has 823 elements
> > > > Group 10 ( Prot-Masses) has 2433 elements
> > > > Group 11 ( Non-Protein) has 108286 elements
> > > > Group 12 ( CL) has 4 elements
> > > > Group 13 ( HOH) has 423 elements
> > > > Group 14 ( SOL) has 107856 elements
> > > > Group 15 ( CA) has 2 elements
> > > > Group 16 ( NA) has 1 elements
> > > > Group 17 ( Other) has 108286 elements
> > > > Select a group: 1
> > > > Selected 1: 'Protein'
> > > > trn version: GMX_trn_file (single precision)
> > > > Reading frame 0 time 0.000 Fatal error: Error: Too many iterations
> > > > in routine JACOBI
> > > >
> > > > How I could resolve the problem???
> > > >
> > > > then I try g_rms:
> > > >
> > > >
> > > > On Friday 20 May 2005 10:42 am, Anthony Cruz wrote:
> > > > > How I could do that??? a new tpr??
> > > > >
> > > > > On Friday 20 May 2005 7:53 am, Xavier Periole wrote:
> > > > > > Anthony Cruz wrote:
> > > > > > >Hi:
> > > > > > >I run a MD of a protein in water. when I try to analyse the trajectory
> > > > > > > with g_rms the program stop by the following error :
> > > > > > >Fatal error: Error: Too many iterations in routine JACOBI
> > > > > > >What could be the cause? How I can resolve the problem???
> > > > > >
> > > > > > That is certainly due to a mismatch between your reference topology and
> > > > > > the content of
> > > > > > the trajectory ... make an topology that fits the trajectory.
> > > > >
> > > > > _______________________________________________
> > > > > gmx-users mailing list
> > > > > gmx-users at gromacs.org
> > > > > http://www.gromacs.org/mailman/listinfo/gmx-users
> > > > > Please don't post (un)subscribe requests to the list. Use the
> > > > > www interface or send it to gmx-users-request at gromacs.org.
> > > >
> > > >
> > > > ------------------------------
> > > >
> > > > Message: 4
> > > > Date: Wed, 25 May 2005 13:48:43 -0500
> > > > From: "Francesco Mercuri" <mercuri at email.com>
> > > > Subject: [gmx-users] Time reversibility and settle
> > > > To: gmx-users at gromacs.org
> > > > Message-ID: <20050525184843.B8250101D9 at ws1-3.us4.outblaze.com>
> > > > Content-Type: text/plain; charset="iso-8859-1"
> > > >
> > > > > > Hello and thank you for your reply.
> > > > > > However, my question about the time reversibility of a leap-frog integrator
> > > > > > with settle (e.g. for TIP3P water molecules) was concerned with the
> > > > > > actual implementation in Gromacs, rather than the analytical point
> > > > > > of view. In other words, is this implementation of leap-frog +
> > > > > > settle still time-reversible in the limit of infinite numerical
> > > > > > accuracy?
> > > > > > This question arises since I tried different numerical
> > > > > > accuracies (single, double and quadruple precision numbers and
> > > > > > operations; for the latter I had to rewrite large parts of the
> > > > > > code...) with, in all cases, similar (and pretty large)
> > > > > > deviations from the "forward" trajectory when reversing the time.
> > > > > > On the other hand, just bypassing the "settle" subroutine
> > > > > > (e.g. by commenting all the code related to the original
> > > > > > implementation of settle by Miyamoto and Kollman), it works as
> > > > > > expected: the algorithm is almost perfectly time-reversible
> > > > > > and the only "noise" is due to numerical inaccuracies, with
> > > > > > deviations comparable, in the three different cases (single,
> > > > > > double, long double) to the accuracy of numbers.
> > > > > > Thus, I'm trying to understand whether just resetting the position
> > > > > > of atoms, as done in the current implementation of "settle",
> > > > > > still leads to a time-reversible MD (in the limit of infinite
> > > > > > machine accuracy).
> > > > > > Any suggestion about that?
> > > > >
> > > > > In the limit is should be time-reversible, unless there is a bug in
> > > > > the settle algorithm, which I consider very unlikely.
> > > > > Have you tried to use shake or lincs (with high precision settings)?
> > > > >
> > > > > There could be some other issues, but when the unconstrained
> > > > > dynamics is reversible you have probably taken care of these:
> > > > > You should not use temperature and pressure coupling (the new
> > > > > version will have reversible Nose-Hoover coupling).
> > > > > When reversing you should take care to use the proper velocities,
> > > > > as x(t), v(t-t/2dt) is stored you need to take the velocities from
> > > > > the next step when reversing.
> > > > >
> > > > > Berk.
> > > >
> > > > Hi again and thanks for your suggestions.
> > > > Still in the search of better time-reversibility properties, I'd like
> > > > to ask a few more questions.
> > > > If I understood correctly, the leap-frog + settle algorithm is implemented
> > > > in Gromacs according to the following steps:
> > > >
> > > > v'(t+dt/2) = v(t-dt/2) + dt*f(t)/m ("uncorrected" velocities)
> > > > r'(t+dt) = r(t) + dt*v'(t+dt/2) (unconstrained move)
> > > > r(t+dt) = settle[r'(t+dt)] (reset of atomic position calling "csettle")
> > > > v(t+dt/2) = [r(t+dt)-r(t)]/dt (corrected velocities)
> > > >
> > > > is that correct?
> > > > In this case, is not very clear to me how is it possible to "reverse" it,
> > > > e.g. restarting from a [r(t),v(t+dt/2)] configuration and inverting the
> > > > sign of dt, since the settle procedure modifies the atomic position in a
> > > > way that is, apparently, not "reversible".
> > > > Is any modification of the velocities needed?
> > > >
> > > > Another question concerns the preparation of restart files containing
> > > > configurations like the one cited above for time-reversing purposes,
> > > > i.e. [r(t),v(t+dt/2)] (instead of the default [r(t),v(t-dt/2)] ).
> > > > I tried the quick and dirty way by just running one more MD step,
> > > > converting the .trr file in ascii format, manually editing the
> > > > ascii file in order to get the [r(t),v(t+dt/2)] configuration,
> > > > pasting into a .gro file and converting back the .gro file into
> > > > a .trr file for restarting. However, even if I modified the routine
> > > > pr_rvecs in txtdump.c in order to have more decimal digits in the
> > > > ascii files, the accuracy seems to be limited to the order of
> > > > 1.e-8, whereas I expected a better accuracy for double precision
> > > > numbers. Is there any better way to do that? Which is the accuracy
> > > > with which numbers are stored in the .trr file (with double
> > > > precision compilation)?
> > > > Thank you again and best regards,
> > > >
> > > > Francesco
> > > >
> > > > --
> > > > ___________________________________________________________
> > > > Sign-up for Ads Free at Mail.com
> > > > http://promo.mail.com/adsfreejump.htm
> > > >
> > > >
> > > >
> > > > ------------------------------
> > > >
> > > > Message: 5
> > > > Date: Wed, 25 May 2005 20:54:30 +0200
> > > > From: "Jordi Camps" <jcamps at lsi.upc.edu>
> > > > Subject: [gmx-users] Adding ions
> > > > To: "'Discussion list for GROMACS users'" <gmx-users at gromacs.org>
> > > > Message-ID: <200505251854.j4PIsNFw028595 at dagon.lsi.upc.edu>
> > > > Content-Type: text/plain; charset="US-ASCII"
> > > >
> > > > Here I come with more problems :-(
> > > >
> > > > I have a protein file with 126 residues. I generated the topology from the
> > > > .gro file with pdb2gmx. Then I configured the box with editconf. Then I
> > > > added the solvent with genbox. At least I wanted to neutralize the system
> > > > with genion, but it gave me this error:
> > > > Fatal error: pbc_dx called before init_pbc
> > > > This time I think that I followed all the standard steps. Do you know which
> > > > could be the cause of this error?
> > > >
> > > > The command was:
> > > > $ genion -s 1agi.water.top -o 1agi.ion.gro -nn 6
> > > > And the last lines before the error are:
> > > > Selected 12: 'SOL'
> > > > Number of (3-atomic) solvent molecules: 7092
> > > > Doing single force calculation...
> > > > Replacing solvent molecule 1626 (atom 6216) with Cl
> > > > Fatal error: pbc_dx called before init_pbc
> > > >
> > > > Thank you in advance,
> > > >
> > > > --
> > > >
> > > > Jordi Camps Puchades
> > > >
> > > > Instituto Nacional de Bioinformatica (INB) Nodo Computacional GNHC-2
> > > > UPC-CIRI
> > > > c/. Jordi Girona 1-3
> > > > Modul C6-E201 Tel. : 934 011 650
> > > > E-08034 Barcelona Fax : 934 017 014
> > > > Catalunya (Spain) e-mail: jcamps at lsi.upc.edu
> > > >
> > > >
> > > >
> > > > ------------------------------
> > > >
> > > > _______________________________________________
> > > > gmx-users mailing list
> > > > gmx-users at gromacs.org
> > > > http://www.gromacs.org/mailman/listinfo/gmx-users
> > > >
> > > >
> > > > End of gmx-users Digest, Vol 13, Issue 66
> > > > *****************************************
> > > >
> > > _______________________________________________
> > > gmx-users mailing list
> > > gmx-users at gromacs.org
> > > http://www.gromacs.org/mailman/listinfo/gmx-users
> > > Please don't post (un)subscribe requests to the list. Use the
> > > www interface or send it to gmx-users-request at gromacs.org.
> > --
> > David.
> > ________________________________________________________________________
> > David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
> > Dept. of Cell and Molecular Biology, Uppsala University.
> > Husargatan 3, Box 596, 75124 Uppsala, Sweden
> > phone: 46 18 471 4205 fax: 46 18 511 755
> > spoel at xray.bmc.uu.se spoel at gromacs.org http://xray.bmc.uu.se/~spoel
> > ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> >
> > _______________________________________________
> > gmx-users mailing list
> > gmx-users at gromacs.org
> > http://www.gromacs.org/mailman/listinfo/gmx-users
> > Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-users-request at gromacs.org.
> >
> _______________________________________________
> gmx-users mailing list
> gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
--
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596, 75124 Uppsala, Sweden
phone: 46 18 471 4205 fax: 46 18 511 755
spoel at xray.bmc.uu.se spoel at gromacs.org http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
More information about the gromacs.org_gmx-users
mailing list