Re: [gmx-users] Had anybody modied the source code of gromacs
On 19/04/2010 3:58 PM, 聂雪川 wrote: Hello gmx-users, I just have one question about modifing the source code of gromacs4.0.7.For I want to control the motion of one certain atom (e.g. atom 200).I added some codes like /if(n==my_atom)xprime[n][d] = x[n][d]+...;/ in the functionstatic /void do_update_md()/ of the file/ gromacs-4.0.7\src\mdlib\update.c./But when I run it Parallel (mpi run),I find that the atom serial number n isn't the real atom serial number n on each CPU.So the code /if(n==my_atom)/ is futile.How can I achieve it. The domain decomposition algorithm defines a mapping from the local state (i.e. subset of topology and atoms) to the global state. You need to harness that, or don't run in parallel. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] CG-MD simulation of protein, always crash with protein
Hi, the rmsd you show is for the protein in vacuum? It goes up to 0.7 nm! That is a very large deformation, is it expected? Isn't your protein freezing and all the kinetic energy going into translational motion. there is a few things you could try. 1- First you should increase the rlilst to at least 1.3 but 1.4 would be better. 2- And/or reduce the nstlist to 5. You should modify at least one of the two, that will increases the cost of your simulation but should make increase energy conservation. Then you should try a smaller time step 0.020 ps might just remove all the lincs warnings. XAvier. On Fri, Apr 16, 2010 at 5:17 PM, Martti Louhivuori m.j.louhivu...@rug.nl wrote: On 15 Apr 2010, at 18:06, Trang wrote: My target system is a protein with lipid molecules added randomly (using GENBOX). Running MD, I expect to I hope you're using a larger van der Waals distance (0.24nm or so) when inserting the lipids. I tried 0.3 - 0.5. It looked fine, but didn't work. broke down the problem, that is, to run md simulation for the protein molecule only, and in vacuum. Still no improvement. Although all the distances in the minimized structure are visually proper, the system exploded. If you can't run the protein even in vacuum, then the problem is either in your MD parameters, starting co-ordinates or some simple mistake somewhere. Since the .mdp files you posted seem ok, my bet is on the starting co-ordinates, just as Xavier proposed. Could you post the system topology (.top) and the protein topology (.itp), so we can rule out any mistakes in those? The protein topo is too large, so I put it here. http://pastie.org/926617 Here is the system topology. #include ../martini_v2.1.itp #include ../martini_v2.1_aminoacids.itp #include 11BHSD1.itp ; Include Position restraint file #ifdef POSRES #include posre.itp #endif [ system ] 11BHSD1 [ molecules ] 11BHSD1 1 You should also double check whether you have close contacts between some atoms in the protein; e.g. in VMD this is easily done using the dynamic bonds representation. -martti- Minimized structure showed no close contact (with distance cutoff 1.6, even to 1.9). Production run stopped at step 175552 (5266.56 ps) with Segmentation fault, preceeded with a lot of Lincs warnings (I set ring_bonds = constraints). This is the mdp file, this file goes along with the system in vaccuum that crashed. I'm not sure if too long simulation time can be the cause. The strange thing is that the system seemed to be relatively stable at the time of crash. rmds.xvg is below, in case you want to see the graph first-hand. -MD--- integrator = md tinit= 0.0 dt = 0.030 nsteps = 90 nstcomm = 1 comm-grps= nstxout = 5000 nstvout = 5000 nstfout = 0 nstlog = 2000 nstenergy= 2000 nstxtcout= 1000 xtc_precision= 100 xtc-grps = energygrps = nstlist = 10 ns_type = grid pbc = xyz rlist= 1.2 coulombtype = Shift rcoulomb_switch = 0.0 rcoulomb = 1.2 epsilon_r= 15 vdw_type = Shift rvdw_switch = 0.9 rvdw = 1.2 DispCorr = No tcoupl = Berendsen tc-grps = Protein tau_t= 0.3 ref_t= 323 Pcoupl = berendsen Pcoupltype = isotropic tau_p= 3.0 compressibility = 3e-5 ref_p= 1.0 gen_vel = no gen_temp = 323 gen_seed = 666 constraints = none constraint_algorithm = Lincs unconstrained_start = no lincs_order = 4 lincs_warnangle = 30 -- -rmsd.xvg--- @title RMSD @xaxis label Time (ps) @yaxis label RMSD (nm) @TYPE xy @ subtitle Protein after lsq fit to Protein 0.0000.0001400 150.0000.6150032 300.0000.6162945 450.0000.6451609 600.0000.6317724 750.0000.6501579 900.0000.6642945 1050.0000.6742513 1200.0000.6651280 1350.0000.6965100 1500.0000.6856629 1650.0000.7271055 1800.0000.7174531 1950.0000.7088170 2100.0000.7492073 2250.0000.7352809 2400.0000.7196461 2550.0000.7162255 2700.0000.7277457 2850.0000.7121301 3000.0000.7084662 3150.0000.7153199 3300.000
Re: [gmx-users] CG-MD simulation of protein, always crash with protein
On 19 Apr 2010, at 06:49, Trang wrote: The protein topo is too large, so I put it here. http://pastie.org/926617 The protein topology is fine, but you could try also with elastic bonds instead of dihedrals in the extended regions (--elastic option). Elastic bonds are more robust than dihedrals, especially if running with a larger timestep (30fs). Here is the system topology. #include ../martini_v2.1.itp #include ../martini_v2.1_aminoacids.itp The last line is unnecessary. Use it only when simulating single aminoacids instead of a polypeptide. Minimized structure showed no close contact (with distance cutoff 1.6, even to 1.9). Production run stopped at That's only 0.16-0.19nm, so overlap is still possible. step 175552 (5266.56 ps) with Segmentation fault, preceeded with a lot of Lincs warnings (I set ring_bonds = constraints). This is the mdp file, this file goes along with the system in vaccuum that crashed. I'm not sure if What's different this time around? Last time you said you couldn't run it even in vacuum. Sounds like you are just hitting statistical limits with your combination of (extended) dihedrals + 30fs timestep. Use a smaller timestep (20fs for example) or switch to local elastic bonds for the beta-strands... and solvate the system (without overlap) so the simulation makes some sense! -martti- -- Post-doctoral research fellow Moleculaire Dynamica University of Groningen Nijenborgh 4, 9747AG Groningen, the Netherlands tel. +(31) 50 363 4339 | fax. +(31) 50 363 4398 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] help
Hi gmx users, I used g_lie command to compute free energy. The command and the output is as follow [shik...@venus 30078-dyn10ns]$ g_lie -f edr.edr -ligand UNK -o lie10ns.xvg -ligand string UNK Name of the ligand in the energy file Opened edr.edr as single precision energy file Using the following energy terms: LJ: LJ-SR:Protein-UNK LJ-LR:Protein-UNK LJ-14:Protein-UNK LJ-SR:SOL-UNK LJ-LR:SOL-UNK LJ-14:SOL-UNK LJ-SR:UNK-rest LJ-LR:UNK-rest LJ-14:UNK-rest Coul: Coul-SR:Protein-UNK Coul-14:Protein-UNK Coul-SR:SOL-UNK Coul-14:SOL-UNK Coul-SR:UNK-rest Coul-14:UNK-rest Last frame read 5 time 1.000 DGbind = -88.393 (5.098) gcq#306: Miggida-Miggida-Miggida-Mac (Kriss Kross) The value is coming DGbind = -88.393 (5.098). Whether this is correct method to get the free energy of binding of a ligand. Any other command or method is there to get the free energy of binding of ligands or can we do the MM/PBSA and MM/GBSA analysis from the gromacs dynamics. -- Shikhar Gupta Senior Research Fellow Pharmacoinformatics Department Block- A (Room No.- 208) National Institute of Pharmaceutical Education Research( NIPER ) Sec- 67, S.A.S Nagar Mohali, Punjab (India) Web-Site: www.niper.ac.in PIN- 160062 Email:shik_...@rediffmail.com,shik...@gmail.com-- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Re: Gromacs issue
Please keep GROMACS correspondence on the mailing list. You stand a much better chance of getting a useful reply that can be saved for the future. On 19/04/2010 10:35 PM, Martin Vartorelli wrote: Hi, I'm a Gromacs user. I have a problem and I couldn't find any answer in the mailing list and I'm very hurry. That's irrelevant to whether you should be pestering an individual who hasn't solicited your correspondence :-) I'm asking you because I think that you already has an issue like mine. No, I've never had an issue like yours. You don't want to irritate the people who might be best place to help you by spamming them. I'm trying to use tabulated potentials for a very simple system: a polymer chain with only 3 kind of atoms: A, B and C, all of it with zero charge. A--(B)n--C The only interactions that accounts are AB, AC, BB and CC, so in my mdp file I have writen: coulombtype = user vdw-type= user energygrps = A B C energygrp_table = AB AC BB BC Also, I'm using tabulated bonded potentials for bonds and angles. The main questions are how to tell to mdrun that he must read my table files and can I be sure that the reading is OK? Make sure you've read the manual and relevant wiki web page thoroughly, and you should experiment on a simple system first. Mark For the bonding part I have the following table files: table_b1.xvg table_a1.xvg table_a2.xvg table_a3.xvg And in the mdrun commandline I'm using the option -tableb table. But how to handle the non bonding table files? table_A_B.xvg table_A_C.xvg table_B_B.xvg table_B_C.xvg Must I write energygrps = A B C energygrp_table = AC BB BC in the mdp file and use the tables table.xvg table_A_C.xvg table_B_B.xvg table_B_C.xvg with the command line for mdrun saying -table table.xvg? Another thing, because the tables doesn't have to contain only zero columns... Can I put anything on the f and f' columns (because the charge of my atoms is zero)? Can I put (for example) 0.5 in the g and g' columns and my potential (and the force) minus 0.5 in columns h and h'? Any help will be appreciated. Thank you. Lic. Martin R. Vartorelli INTEC-UNL-CONICET (Institute of Technological Development for the Chemical Industry) Güemes 3450 3000 Santa Fe, Argentina Phone: +54 (342) 455-9174 (ext. 2072, office PPB23) -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Problem with tabulated potentials for 3 different atoms
I'm trying to use tabulated potentials for a very simple system: a polymer chain with only 3 kind of atoms: A, B and C, all of it with zero charge. A--(B)n--C The only interactions that accounts are AB, AC, BB and CC, so in my mdp file I have writen: coulombtype = user vdw-type = user energygrps = A B C energygrp_table = AB AC BB BC Also, I'm using tabulated bonded potentials for bonds and angles. The main questions are how to tell to mdrun that he must read my table files and can I be sure that the reading is OK? For the bonding part I have the following table files: table_b1.xvg table_a1.xvg table_a2.xvg table_a3.xvg And in the mdrun commandline I'm using the option -tableb table. But how to handle the non bonding table files? table_A_B.xvg table_A_C.xvg table_B_B.xvg table_B_C.xvg Must I write energygrps = A B C energygrp_table = AC BB BC in the mdp file and use the tables table.xvg table_A_C.xvg table_B_B.xvg table_B_C.xvg with the command line for mdrun saying -table table.xvg? Another thing, because the tables doesn't have to contain only zero columns... Can I put anything on the f and f' columns (because the charge of my atoms is zero)? Can I put (for example) 0.5 in the g and g' columns and my potential (and the force) minus 0.5 in columns h and h'? I'm using the 4.0.7 version. Again as an abstract: System: 1 chain A--(B)n--C Bonds: A-B = B-B = B-C -- only one table: table_b1.xvg Angles: A-B-B, B-B-B, B-B-C -- 3 tables: table_a1.xvg, table_a2.xvg, table_a3.xvg Charges: none VDW interactions: A-B, A-C, B-B, B-C -- 4 tables table_A_B.xvg, table_A_C.xvg, table_B_B.xvg, table_B_C.xvg *.mdp lines: coulombtype = user vdw-type = user energygrps = A B C energygrp_table = ? run commandline: mdrun -s -c -e -x -g -table ? -tableb ? Any help will be appreciated. Thank you. Martin R. Vartorelli -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] about H-Bonds g_hbond
Hi all I am getting an error in g_hbond such as Program g_hbond_d, VERSION 4.0.7 Source code file: gmx_hbond.c, line: 565 Fatal error: Error in func_type Position Rest. Could anyone tell me how to solve this problem. thanks in advance E R Azhagiya singam Send free SMS to your Friends on Mobile from your Yahoo! Messenger. Download Now! http://messenger.yahoo.com/download.php-- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] about H-Bonds g_hbond
babu gokul wrote: Hi all I am getting an error in g_hbond such as Program g_hbond_d, VERSION 4.0.7 Source code file: gmx_hbond.c, line: 565 Fatal error: Error in func_type Position Rest. Could anyone tell me how to solve this problem. thanks in advance If you have position restraints defined, you can't run g_hbond for some reason. The work-around is to re-generate a .tpr file that does not use position restraints and do the analysis. -Justin E R Azhagiya singam Send free SMS to your Friends on Mobile from your Yahoo! Messenger. Download Now! http://messenger.yahoo.com/download.php -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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/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/mailing_lists/users.php
Re: [gmx-users] Problem with tabulated potentials for 3 different atoms
Dear Martin As Mark has already told you - if you have a problem with gromacs email the list and wait for a reply. Don't just email people you don't know from Adam and expect them to reply because you are in a rush. The file that you have read from the mailing list contains everything I know about tabulated potentials in gromacs. I can't say I fully understand your email but I think that your questions (other than the stuff on bonded potentials) are all covered on the wikki page and certainly are covered in the pdf file I wrote. In short For the bonding part I have the following table files: table_b1.xvg table_a1.xvg table_a2.xvg table_a3.xvg I know nothing about using tabulated potentials for bonded interactions. If this is what is in the manual though its probably right. Must I write energygrps = A B C energygrp_table = AC BB BC in the mdp file and use the tables table.xvgThis has the A B interaction in? table_A_C.xvg table_B_B.xvg table_B_C.xvg with the command line for mdrun saying -table table.xvg? Yes this is exactly what the wikki tells you to do. Why do you think it wont work in your case? Incidentally, you can work out if this is working by doing one md step in gromacs and comparing the energy you obtain with the energy you get from another code, which you are confident produces correct energies. Another thing, because the tables doesn't have to contain only zero columns... Can I put anything on the f and f' columns (because the charge of my atoms is zero)? Yes. I always use f=1/r even though I know it will be ignored (with this choice f'=1/r**2 obviously) Can I put (for example) 0.5 in the g and g' columns and my potential (and the force) minus 0.5 in columns h and h'? I don't understand what you want to do here. If you put 0.5 in g then g'=0, which you can work out using, what I assume, is the first thing you learnt about calculus. ciao Gareth -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Directional calculation Sg and Potential
Hello, This is the second time I am sending this message. I have a protein which has a few Thr residues on one side of itself. I would like to know how should I setup mt index file and use which GROMACS tools to calculate potential energy and Sg of water molecules that only face these residues as a function of space coordinate such as z (not radial coordinate r), assuming I want my Thr residues be at the z=0 of the profile of Sg vs. z or PE vs z. Thanks for your helps. Paymon -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Re: manual eq. 4.74-4.75 (dihedral restraints) head scratcher
Justin,
I appreciate your helpful reply. To summarize the additional points I
make below:
(1) the code is right as far as calculating the dihedral restraint
potentials, but the manual has a bug, and
(2) the settings for dihedral restraints should be documented,
preferably in the wiki, which I'm not able to edit.
Justin wrote:
Daniel L. Ensign wrote:
Hello gmx-users, you rock and rollers,
Equations 4.74 and 4.75 in my copy of the manual have (please pardon my
pseudo-LaTeX):
(4.74)
\phi' = (\phi - \phi_0) MOD 2\pi
(4.75)
V(\phi') = \frac{1}{2}k ( \phi' - \phi_0 - \Delta \phi )^2 if \phi'
\Delta \phi
or
V(\phi') = 0 if \phi' \leq \Delta \phi
but should there be absolute values around all of the \phi-\phi_0? Which
way is it in the code -- with the absolute distance between \phi and
\phi_0 or the directed distance?
It looks like absolute values are considered. From src/gmxlib/dihres.c:
/* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
* force changes if we just apply a normal harmonic.
* Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
* This means we will never have the periodicity problem, unless
* the dihedral is Pi away from phiO, which is very unlikely due to
* the potential.
*/
dp = phi-phi0;
if (fabs(dp) dphi) {
/* dp cannot be outside (-2*pi,2*pi) */
if (dp = M_PI)
dp -= 2*M_PI;
else if(dp -M_PI)
dp += 2*M_PI;
Thanks for pointing that out. As I look closer, I also see
if (dp dphi)
ddp = dp-dphi;
else if (dp -dphi)
ddp = dp+dphi;
else
ddp = 0;
followed by
vtot += 0.5*kfac*ddp*ddp;
Your code snippet indicates that the code is correct (absolute value
of dp when calculating the angle modulus) and mine also indicates that
the code is correct (absolute value of dp when calculating the
restraint energy).
The manual, however, should have absolute values in each place that
phi-phi0 appears. How can the manual be corrected expending the least
effort?
Also, as far as I can tell (and some mornings I definitely don't read
too good) neither the manual nor
http://wiki.gromacs.org/index.php/Dihedral_Restraints define the fields
in [ dihedral_restraints ], although the latter does name them. There, I
see
[ dihedral_restraints ]
; ai ajakal type label phi dphi kfac power
57 915 1 1 180 0 1 2
ai, aj, ak, al = atom numbers, obviously
type = ?, but I'm guessing there's only one type anyway
Probably so.
label = what is this one?
Looks to be bookkeeping. The code doesn't seem to use it other than to print
debug information, but I could be wrong since I haven't surfed
around it very long.
snip
kfac is the force constant, probably
Indirectly. This term is equivalent to the fac value in distance restraints.
Since the force constant is specified in the .mdp file, different restraints
would otherwise have to be restrained with equivalent force constants. The
value of kfac is multiplied by the value of dihre_fc in the .mdp
file, so that
different restraints could have different force constants.
So then it seems there are two ways to set force constants -- by
setting dihre_fc in the mdp file and by setting kfac in the dihedral
restraints itp file. Good to know.
power = what is this one? Does 2 give me harmonic constraints?
Not a clue on this one. Also doesn't seem to be used in the code, but maybe
it's somewhere outside of dihres.c.
For now I'll just cross my fingers and put a 2 ... but it would be
good for all of these doohickeys to be documented.
I'd be happy to translate any answers given to the wiki, assuming I get
answered and that I'm allowed to edit the wiki.
sniip
I might get around to updating this page with the above information, unless
there is more information to be considered, or I'm wrong :)
I don't appear to have the appropriate privileges for editing, but
that's okay, as it seems like a lot of responsibility.
Have fun,
Dan
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Re: manual eq. 4.74-4.75 (dihedral restraints) head scratcher
Daniel L. Ensign wrote:
Justin,
I appreciate your helpful reply. To summarize the additional points I
make below:
(1) the code is right as far as calculating the dihedral restraint
potentials, but the manual has a bug, and
Yes, perhaps some additional information in the manual would be useful.
(2) the settings for dihedral restraints should be documented,
preferably in the wiki, which I'm not able to edit.
I've updated it:
http://www.gromacs.org/Documentation/How-tos/Dihedral_Restraints
-Justin
Justin wrote:
Daniel L. Ensign wrote:
Hello gmx-users, you rock and rollers,
Equations 4.74 and 4.75 in my copy of the manual have (please pardon my
pseudo-LaTeX):
(4.74)
\phi' = (\phi - \phi_0) MOD 2\pi
(4.75)
V(\phi') = \frac{1}{2}k ( \phi' - \phi_0 - \Delta \phi )^2 if \phi'
\Delta \phi
or
V(\phi') = 0 if \phi' \leq \Delta \phi
but should there be absolute values around all of the \phi-\phi_0? Which
way is it in the code -- with the absolute distance between \phi and
\phi_0 or the directed distance?
It looks like absolute values are considered. From src/gmxlib/dihres.c:
/* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
* force changes if we just apply a normal harmonic.
* Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
* This means we will never have the periodicity problem, unless
* the dihedral is Pi away from phiO, which is very unlikely due to
* the potential.
*/
dp = phi-phi0;
if (fabs(dp) dphi) {
/* dp cannot be outside (-2*pi,2*pi) */
if (dp = M_PI)
dp -= 2*M_PI;
else if(dp -M_PI)
dp += 2*M_PI;
Thanks for pointing that out. As I look closer, I also see
if (dp dphi)
ddp = dp-dphi;
else if (dp -dphi)
ddp = dp+dphi;
else
ddp = 0;
followed by
vtot += 0.5*kfac*ddp*ddp;
Your code snippet indicates that the code is correct (absolute value of
dp when calculating the angle modulus) and mine also indicates that the
code is correct (absolute value of dp when calculating the restraint
energy).
The manual, however, should have absolute values in each place that
phi-phi0 appears. How can the manual be corrected expending the least
effort?
Also, as far as I can tell (and some mornings I definitely don't read
too good) neither the manual nor
http://wiki.gromacs.org/index.php/Dihedral_Restraints define the fields
in [ dihedral_restraints ], although the latter does name them. There, I
see
[ dihedral_restraints ]
; ai ajakal type label phi dphi kfac power
57 915 1 1 180 0 1 2
ai, aj, ak, al = atom numbers, obviously
type = ?, but I'm guessing there's only one type anyway
Probably so.
label = what is this one?
Looks to be bookkeeping. The code doesn't seem to use it other than
to print
debug information, but I could be wrong since I haven't surfed around
it very long.
snip
kfac is the force constant, probably
Indirectly. This term is equivalent to the fac value in distance
restraints.
Since the force constant is specified in the .mdp file, different
restraints
would otherwise have to be restrained with equivalent force
constants. The
value of kfac is multiplied by the value of dihre_fc in the .mdp
file, so that
different restraints could have different force constants.
So then it seems there are two ways to set force constants -- by setting
dihre_fc in the mdp file and by setting kfac in the dihedral restraints
itp file. Good to know.
power = what is this one? Does 2 give me harmonic constraints?
Not a clue on this one. Also doesn't seem to be used in the code, but
maybe
it's somewhere outside of dihres.c.
For now I'll just cross my fingers and put a 2 ... but it would be good
for all of these doohickeys to be documented.
I'd be happy to translate any answers given to the wiki, assuming I get
answered and that I'm allowed to edit the wiki.
sniip
I might get around to updating this page with the above information,
unless
there is more information to be considered, or I'm wrong :)
I don't appear to have the appropriate privileges for editing, but
that's okay, as it seems like a lot of responsibility.
Have fun,
Dan
--
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
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/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/mailing_lists/users.php
[gmx-users] Re : keeping the protein in the centre of the simulation box.
Hi all , I have performed a 3ns simulation and after that I have found that the protein has moved out of the simulation box. Can anybody tell how can I fix this problem without performing the simulation again .. Thanks -- Bharat M.Sc. Bioinformatics (Final year) Centre for Bioinformatics Pondicherry University Puducherry India Mob. +919962670525 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Re : keeping the protein in the centre of the simulation box.
bharat gupta wrote: Hi all , I have performed a 3ns simulation and after that I have found that the protein has moved out of the simulation box. Can anybody tell how can I fix this problem without performing the simulation again .. Thanks One of the most commonly-asked non-problems across this list... http://www.gromacs.org/Documentation/Terminology/Periodic_Boundary_Conditions -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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/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/mailing_lists/users.php
[gmx-users] reference for gromos 87 within PRODRG
Hi, I have generated some topology files using the gromos87 forcefield within the PRODRG server. Does anyone have a reference for this force-field? I can't seem to find it, Thanks Jenny -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] reference for gromos 87 within PRODRG
Jennifer Williams wrote: Hi, I have generated some topology files using the gromos87 forcefield within the PRODRG server. Does anyone have a reference for this force-field? I can't seem to find it, There is some useful information here at the bottom of the page, as well as the link to the archive post regarding GROMOS87/ffgmx. I think the original reference is based out of the GROMOS87 user's manual, which is not in the public domain. http://www.gromacs.org/Documentation/Terminology/Force_Fields/GROMOS -Justin Thanks Jenny -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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/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/mailing_lists/users.php
Re: [gmx-users] Directional calculation Sg and Potential
On 20/04/2010 1:16 AM, Paymon Pirzadeh wrote: Hello, This is the second time I am sending this message. That'd be because nobody had a good answer the first time... I have a protein which has a few Thr residues on one side of itself. I would like to know how should I setup mt index file and use which GROMACS tools to calculate potential energy and Sg of water molecules that only face these residues as a function of space coordinate such as z (not radial coordinate r), assuming I want my Thr residues be at the z=0 of the profile of Sg vs. z or PE vs z. This sounds like a complex analysis that you will need to be prepared to do a lot of work to code yourself. Cunning use of index groups with (for example) trjconv, g_traj, trjorder, mdrun -rerun and whatever tool calculates Sg will be required. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Problem with tabulated potentials
Hello Gromacs users, This is constructive criticism to some people on the list... Tabulated potentials are an advanced feature of Gromacs and the people who is looking for help on the mailing list (certainly not experienced) deserve a little respect, simply because they eventually will contribute to the human knowledge, like probably every single person that search answers on the list. I'm saying this because of the dramatic answers that can be readed here: http://lists.gromacs.org/pipermail/gmx-users/2010-April/050256.html Gromacs manual requieres: knowledge about molecular modeling in general and patience. More than a guide the manual is a challenge and specially the tabulated potentials stuff. People often doesn't have time to learn a new programming languaje to see inside the code to understand how the machine works (this is the reason for a manual!). Because I usually work with high MW polymers, I know the feeling of discouragement when you're trying to use Gromacs for this kind of systems. But If I had received such an hysterical answer like http://lists.gromacs.org/pipermail/gmx-users/2010-April/050256.html or http://lists.gromacs.org/pipermail/gmx-users/2010-April/050252.html without any doubt, I'd change the simulation software. Mark.Abraham: take an anxiolytic pill if you continue suffering this mood disorder when you receive an e-mail with a simple consultation. Gareth Tribello: your guide doesn't tell anything that complements the manual, why you have written it? (go to the talk shows if you want to start to be famous) I had to eat my pasta cold because of this gossips! Peace. J.M.S. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Problem with tabulated potentials
On 20/04/2010 7:32 AM, Jhony Tolengo wrote: Hello Gromacs users, This is constructive criticism to some people on the list... Suggesting that people have delusions of grandeur or need medication is not constructive criticism. That's abuse, and has no place on this list. Tabulated potentials are an advanced feature of Gromacs and the people who is looking for help on the mailing list (certainly not experienced) deserve a little respect, simply because they eventually will contribute to the human knowledge, like probably every single person that search answers on the list. They get respect, but they will get more respect (and more help!) if they participate actively in their learning process. The advice here http://www.catb.org/%7Eesr/faqs/smart-questions.html is excellent reading for people participating on mailing lists. I'm saying this because of the dramatic answers that can be readed here: http://lists.gromacs.org/pipermail/gmx-users/2010-April/050256.html There's nothing dramatic or hysterical there. Gromacs manual requieres: knowledge about molecular modeling in general and patience. More than a guide the manual is a challenge and specially the tabulated potentials stuff. People often doesn't have time to learn a new programming languaje to see inside the code to understand how the machine works (this is the reason for a manual!). Because I usually work with high MW polymers, I know the feeling of discouragement when you're trying to use Gromacs for this kind of systems. Doing good science is difficult and time-consuming and MD is typical here. Using advanced features of complex software is always time-consuming. GROMACS is excellent in that it makes it possible to use tabulated potentials without needing to write or understand code. Examples of such tables are provided in the code base. If you really wish to give constructive criticism, what does the manual lack in regard to usage of tables? How should Gareth and Justin's web contribution be improved? Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Problem with tabulated potentials
The only hysterical email I've seen is the one I just read. Posting etiquette tips to users is not unreasonable. For those of us who receive dozens of personal requests a week, it becomes somewhat tiresome. For this reason, we must use very clear language that directs the user to the proper resources. I have experienced several communications where users simply do not understand (for a number of reasons, usually a language barrier or stubbornness) that I am unable to help them until I make comments like Mark's. I would ask you, what have you contributed? If you can improve upon any resources, please do. Gromacs is a user-supported community. If you think someone is wrong or you can explain something better, step up. Otherwise, don't undermine someone else's honest attempt to help out. I think Gareth's document is a nice walk-through of a simple example of tabulated potentials, which complements the three or so paragraphs in the manual. Mark has helped literally thousands of people over this list, and his advice (as well as several others') was invaluable to me in preparing a difficult project. If you have no other comments than to bring others down, feel free to delete the emails you don't care about. This list is for productive discussion and scholarly exchange. No need for these snide remarks. None of what you've posted is constructive, as you claim, but your mail is now archived as a record of an unprofessional and unnecessary rant. -Justin Jhony Tolengo wrote: Hello Gromacs users, This is constructive criticism to some people on the list... Tabulated potentials are an advanced feature of Gromacs and the people who is looking for help on the mailing list (certainly not experienced) deserve a little respect, simply because they eventually will contribute to the human knowledge, like probably every single person that search answers on the list. I'm saying this because of the dramatic answers that can be readed here: http://lists.gromacs.org/pipermail/gmx-users/2010-April/050256.html Gromacs manual requieres: knowledge about molecular modeling in general and patience. More than a guide the manual is a challenge and specially the tabulated potentials stuff. People often doesn't have time to learn a new programming languaje to see inside the code to understand how the machine works (this is the reason for a manual!). Because I usually work with high MW polymers, I know the feeling of discouragement when you're trying to use Gromacs for this kind of systems. But If I had received such an hysterical answer like http://lists.gromacs.org/pipermail/gmx-users/2010-April/050256.html or http://lists.gromacs.org/pipermail/gmx-users/2010-April/050252.html without any doubt, I'd change the simulation software. Mark.Abraham: take an anxiolytic pill if you continue suffering this mood disorder when you receive an e-mail with a simple consultation. Gareth Tribello: your guide doesn't tell anything that complements the manual, why you have written it? (go to the talk shows if you want to start to be famous) I had to eat my pasta cold because of this gossips! Peace. J.M.S. -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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/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/mailing_lists/users.php
Re: [gmx-users] Concerns with g_wham
Hello again, As you mentioned, it makes sense that the starting umbrella potential a in the equation U = (1/2)K(x-a)^2 should be given in the .tpr file. I have looked over the manual and online though, and I still cannot see what portion of the pull code allows me to define the a. It seems that maybe GROMACS just takes the starting distance between the two atoms, and centers the umbrella potential around that distance? For instance, my starting configuration starts with the reference group and pull group 0.28 nm away. I have turned on the pull code in the .tpr file in order to do an umbrella sample and it looks like this: pull = umbrella pull_geometry = distance pull_dim = Y Y Y pull_start = yes pull_ngroups = 1 pull_group0 = Na+ pull_group1 = I- pull_rate1 = 0.00 pull_k1 = 5000 I don't see how any of these designate where the bottom of the umbrella potential well lies (except for maybe pull_start - perhaps setting it to yes puts centers the umbrella potential on the COM of the system). So instead, maybe GROMACS realizes that Na+ and I- start 0.28 nm away, and uses this value as the center of the umbrella potential. Is this the case? Or do one of these pull code designations dictate a? Thanks, Jennifer On Fri, Apr 16, 2010 at 4:20 AM, Jochen Hub joc...@xray.bmc.uu.se wrote: Jennifer Casey wrote: Hello, I have been using g_wham, but I have a few questions that I can't find answers to online. When using WHAM, one does not need the forces between the pull groups to calculate the PMF, yet g_wham won't run without it. Is there a reason for this? Hi Jennifer! As pointed out by Chris, you can use provide the pull positions (g_wham -ix) OR the pull forces with -if. When giving the pullf files, g_wham simply computes the pull positions from the forces and calculates the histograms etc... Also, when using the pull code, I am allowed to define the spring constant K for umbrella sampling, but I do not designate where the umbrella potential is centered. How does gromacs determine this? Umbrella positions and force constants (and pull geometry) are taken from the tpr files. That's why they must be provided to g_wham. I am interested as I would like to create a PMF using umbrella integration (from code I will write myself) rather than use WHAM. To do this and still use the umbrella sampling runs used with GROMACS, I need to know where my umbrella potentials are centered. If you want to integrate the mean forces, I would not use pull=umbrella but pull=constraint. You should get the same PMF compared to umbrella sampling and using g_wham. Please let me know if g_wham makes any trouble. Cheers, Jochen Thank you, Jennifer -- --- Dr. Jochen Hub Molecular Biophysics group Dept. of Cell Molecular Biology Uppsala University. Box 596, 75124 Uppsala, Sweden. Phone: +46-18-4714451 Fax: +46-18-511755 --- -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Concerns with g_wham
Jennifer Casey wrote: Hello again, As you mentioned, it makes sense that the starting umbrella potential a in the equation U = (1/2)K(x-a)^2 should be given in the .tpr file. I have looked over the manual and online though, and I still cannot see what portion of the pull code allows me to define the a. It seems that maybe GROMACS just takes the starting distance between the two atoms, and centers the umbrella potential around that distance? For instance, my starting configuration starts with the reference group and pull group 0.28 nm away. I have turned on the pull code in the .tpr file in order to do an umbrella sample and it looks like this: pull = umbrella pull_geometry = distance pull_dim = Y Y Y pull_start = yes pull_ngroups = 1 pull_group0 = Na+ pull_group1 = I- pull_rate1 = 0.00 pull_k1 = 5000 I don't see how any of these designate where the bottom of the umbrella potential well lies (except for maybe pull_start - perhaps setting it to yes puts centers the umbrella potential on the COM of the system). So instead, maybe GROMACS realizes that Na+ and I- start 0.28 nm away, and uses this value as the center of the umbrella potential. Is this the case? Or do one of these pull code designations dictate a? You've got it right. The parameter pull_start dictates whether or not the starting distance is added to the value of pull_init1. The default for pull_init1 is zero, so setting pull_start = yes tells grompp (and then mdrun) to add the COM distance of the starting conformation to pull_init (per the manual). So the position of the umbrella is 0.28 (pull_start) + 0 (pull_init1) = 0.28 nm, in your case. -Justin Thanks, Jennifer On Fri, Apr 16, 2010 at 4:20 AM, Jochen Hub joc...@xray.bmc.uu.se mailto:joc...@xray.bmc.uu.se wrote: Jennifer Casey wrote: Hello, I have been using g_wham, but I have a few questions that I can't find answers to online. When using WHAM, one does not need the forces between the pull groups to calculate the PMF, yet g_wham won't run without it. Is there a reason for this? Hi Jennifer! As pointed out by Chris, you can use provide the pull positions (g_wham -ix) OR the pull forces with -if. When giving the pullf files, g_wham simply computes the pull positions from the forces and calculates the histograms etc... Also, when using the pull code, I am allowed to define the spring constant K for umbrella sampling, but I do not designate where the umbrella potential is centered. How does gromacs determine this? Umbrella positions and force constants (and pull geometry) are taken from the tpr files. That's why they must be provided to g_wham. I am interested as I would like to create a PMF using umbrella integration (from code I will write myself) rather than use WHAM. To do this and still use the umbrella sampling runs used with GROMACS, I need to know where my umbrella potentials are centered. If you want to integrate the mean forces, I would not use pull=umbrella but pull=constraint. You should get the same PMF compared to umbrella sampling and using g_wham. Please let me know if g_wham makes any trouble. Cheers, Jochen Thank you, Jennifer -- --- Dr. Jochen Hub Molecular Biophysics group Dept. of Cell Molecular Biology Uppsala University. Box 596, 75124 Uppsala, Sweden. Phone: +46-18-4714451 Fax: +46-18-511755 --- -- gmx-users mailing listgmx-users@gromacs.org mailto:gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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 mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee 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/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/mailing_lists/users.php
[gmx-users] Histrionic answers doesn't help!
Mark Abraham and Gareth Tribello: Because of your embarrassing answers (this is abuse!) to mailing list users, many of them simply will choose other sotware than Gromacs to simulate systems. For you health, please don't reply, take both an anxiolytic of any kind or take vacations! Peace. J.M.S. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] help
Respected Sir/Madam, I am student of Tribhuwan University and trying to find out diffusion coefficient of D2O in H2O. But due to lack of access of journal i am not being able to complete my job. Can any tell me how to define angle in the topology file and the value of harmonic constant of D2O. I am really in very serious trouble. Udaya Raj Dahal Tribhuwan University Nepal _ Your E-mail and More On-the-Go. Get Windows Live Hotmail Free. https://signup.live.com/signup.aspx?id=60969-- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
RE: [gmx-users] help
See section 5.7.1 of the manual , it shows you exactly how you define angles, bond, dihedrals etc for any type of molecule in the topology files for GROMACS. If after the theory, see section 4.2 Catch ya, Dr. Dallas Warren Drug Delivery, Disposition and Dynamics Monash Institute of Pharmaceutical Sciences, Monash University 381 Royal Parade, Parkville VIC 3010 dallas.war...@pharm.monash.edu.au +61 3 9903 9167 - When the only tool you own is a hammer, every problem begins to resemble a nail. From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] On Behalf Of udaya dahal Sent: Tuesday, 20 April 2010 10:47 AM To: gmx-users@gromacs.org Subject: [gmx-users] help Respected Sir/Madam, I am student of Tribhuwan University and trying to find out diffusion coefficient of D2O in H2O. But due to lack of access of journal i am not being able to complete my job. Can any tell me how to define angle in the topology file and the value of harmonic constant of D2O. I am really in very serious trouble. Udaya Raj Dahal Tribhuwan University Nepal Your E-mail and More On-the-Go. Get Windows Live Hotmail Free. Sign up now. https://signup.live.com/signup.aspx?id=60969 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Force output
Dear GROMACS developers and users, I have a question regarding GROMACS force output. I need the total force on each atom at every timestep in a constant temperature MD trajectory. I am using the Nose-Hoover thermostat for constant temperature simulations. I know that by using the 'nstfout' keyword in the .mdp file I can tell GROMACS to output forces in the trajectory (.trr) file. Now my question is what exactly is written in the output file as force when I use the above mentioned option. I understand that at first this question may seem a little strange and unnecessary, since everybody knows what a force is, but in the context of Nose-Hoover dynamics there is a slight chance of potential confusion as discussed below. The force information is very crucial for our subsequent analysis, that is why we wanted to make sure that we know exactly what is being printed when one uses the 'nstfout' keyword. Total force on a particle can be interpreted as either the negative gradient of the total potential energy with respect to the position of the atom in question OR the instantaneous rate of change of momentum of the atom. Now for normal Hamiltonian dynamics both of these definitions will give the same number but for the Nose-Hoover dynamics they will be different. Because in the Nose-Hoover equations the rate of change of momentum of an atom is a sum of two terms, first term is the negative gradient of total potential and the second term is momentum of the particle multiplied by a thermostat parameter. So, in Nose-Hoover equation instantaneous rate of change of momentum is NOT equal to the negative gradient of the potential with respect to the atom position. In the light of above discussion can anyone please tell me when GROMACS prints forces in the trajectory file during a Nose-Hoover dynamics, which quantity does it print? Is it the negative gradient of the total potential with respect to the atom position OR the total instantaneous rate of change of momentum which includes an extra term? Thanks a lot, Avisek -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Histrionic answers doesn't help!
Unfortunately Jhony, many users want a black box solution that magically makes answers appear without any grind. Few members in this community have the kind of time required to provide these solutions. Mark's reply was somewhat justified by the nature of the question. ~j Quoting Jhony Tolengo polymersp...@gmail.com: Mark Abraham and Gareth Tribello: Because of your embarrassing answers (this is abuse!) to mailing list users, many of them simply will choose other sotware than Gromacs to simulate systems. For you health, please don't reply, take both an anxiolytic of any kind or take vacations! Peace. J.M.S. -- Jared James Thompson Department of Medicinal Chemistry and Molecular Pharmacology Laboratory for Computational Drug Design and Biology RHPH 504C Heine Pharmacy Building 575 Stadium Mall Drive West Lafayette, IN 47907-2091 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] CG-MD simulation of protein, always crash with protein
On Mon, Apr 19, 2010 at 4:08 PM, XAvier Periole x.peri...@rug.nl wrote: Hi, the rmsd you show is for the protein in vacuum? It goes up to 0.7 nm! That is a very large deformation, is it expected? Isn't your protein freezing and all the kinetic energy going into translational motion. there is a few things you could try. 1- First you should increase the rlilst to at least 1.3 but 1.4 would be better. 2- And/or reduce the nstlist to 5. You should modify at least one of the two, that will increases the cost of your simulation but should make increase energy conservation. Then you should try a smaller time step 0.020 ps might just remove all the lincs warnings. XAvier. I made all the modifications suggested, but the system still crashed (with segmentation fault). However, 'elastic bonds' representation that Martti suggested works and I may need these modifications later. Thanks for being so patient with my problem. On Fri, Apr 16, 2010 at 5:17 PM, Martti Louhivuori m.j.louhivu...@rug.nlwrote: On 15 Apr 2010, at 18:06, Trang wrote: My target system is a protein with lipid molecules added randomly (using GENBOX). Running MD, I expect to I hope you're using a larger van der Waals distance (0.24nm or so) when inserting the lipids. I tried 0.3 - 0.5. It looked fine, but didn't work. broke down the problem, that is, to run md simulation for the protein molecule only, and in vacuum. Still no improvement. Although all the distances in the minimized structure are visually proper, the system exploded. If you can't run the protein even in vacuum, then the problem is either in your MD parameters, starting co-ordinates or some simple mistake somewhere. Since the .mdp files you posted seem ok, my bet is on the starting co-ordinates, just as Xavier proposed. Could you post the system topology (.top) and the protein topology (.itp), so we can rule out any mistakes in those? The protein topo is too large, so I put it here. http://pastie.org/926617 Here is the system topology. #include ../martini_v2.1.itp #include ../martini_v2.1_aminoacids.itp #include 11BHSD1.itp ; Include Position restraint file #ifdef POSRES #include posre.itp #endif [ system ] 11BHSD1 [ molecules ] 11BHSD1 1 You should also double check whether you have close contacts between some atoms in the protein; e.g. in VMD this is easily done using the dynamic bonds representation. -martti- Minimized structure showed no close contact (with distance cutoff 1.6, even to 1.9). Production run stopped at step 175552 (5266.56 ps) with Segmentation fault, preceeded with a lot of Lincs warnings (I set ring_bonds = constraints). This is the mdp file, this file goes along with the system in vaccuum that crashed. I'm not sure if too long simulation time can be the cause. The strange thing is that the system seemed to be relatively stable at the time of crash. rmds.xvg is below, in case you want to see the graph first-hand. -MD--- integrator = md tinit= 0.0 dt = 0.030 nsteps = 90 nstcomm = 1 comm-grps= nstxout = 5000 nstvout = 5000 nstfout = 0 nstlog = 2000 nstenergy= 2000 nstxtcout= 1000 xtc_precision= 100 xtc-grps = energygrps = nstlist = 10 ns_type = grid pbc = xyz rlist= 1.2 coulombtype = Shift rcoulomb_switch = 0.0 rcoulomb = 1.2 epsilon_r= 15 vdw_type = Shift rvdw_switch = 0.9 rvdw = 1.2 DispCorr = No tcoupl = Berendsen tc-grps = Protein tau_t= 0.3 ref_t= 323 Pcoupl = berendsen Pcoupltype = isotropic tau_p= 3.0 compressibility = 3e-5 ref_p= 1.0 gen_vel = no gen_temp = 323 gen_seed = 666 constraints = none constraint_algorithm = Lincs unconstrained_start = no lincs_order = 4 lincs_warnangle = 30 -- -rmsd.xvg--- @title RMSD @xaxis label Time (ps) @yaxis label RMSD (nm) @TYPE xy @ subtitle Protein after lsq fit to Protein 0.0000.0001400 150.0000.6150032 300.0000.6162945 450.0000.6451609 600.0000.6317724 750.0000.6501579 900.0000.6642945
Re: [gmx-users] CG-MD simulation of protein, always crash with protein
Elastic bonds did help the system survive. I'm not really understand why do we need more robust bonds in this system. Does it mean that there exists large forces that cause my system to blow up? Would you mind giving an explanation? Thanks indeed for your taking time to help. Trang On Mon, Apr 19, 2010 at 4:26 PM, Martti Louhivuori m.j.louhivu...@rug.nlwrote: On 19 Apr 2010, at 06:49, Trang wrote: The protein topo is too large, so I put it here. http://pastie.org/926617 The protein topology is fine, but you could try also with elastic bonds instead of dihedrals in the extended regions (--elastic option). Elastic bonds are more robust than dihedrals, especially if running with a larger timestep (30fs). Here is the system topology. #include ../martini_v2.1.itp #include ../martini_v2.1_aminoacids.itp The last line is unnecessary. Use it only when simulating single aminoacids instead of a polypeptide. Minimized structure showed no close contact (with distance cutoff 1.6, even to 1.9). Production run stopped at That's only 0.16-0.19nm, so overlap is still possible. step 175552 (5266.56 ps) with Segmentation fault, preceeded with a lot of Lincs warnings (I set ring_bonds = constraints). This is the mdp file, this file goes along with the system in vaccuum that crashed. I'm not sure if What's different this time around? Last time you said you couldn't run it even in vacuum. Sounds like you are just hitting statistical limits with your combination of (extended) dihedrals + 30fs timestep. Use a smaller timestep (20fs for example) or switch to local elastic bonds for the beta-strands... and solvate the system (without overlap) so the simulation makes some sense! -martti- -- Post-doctoral research fellow Moleculaire Dynamica University of Groningen Nijenborgh 4, 9747AG Groningen, the Netherlands tel. +(31) 50 363 4339 | fax. +(31) 50 363 4398 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
