Re: [gmx-users] Normal-mode analysis in Gromacs
Thank you David for your reply. In fact, I have already been using the double precision version of Gromacs. I will greatly appreciate it if other questions could also be covered. Thanks, Kevin Den 2018-03-08 kl. 00:53, skrev Kevin C Chan: > >* Dear Users, *> >* I have encountered several questions when I try to > perform NMA on a protein *>* structure from an equilibrated MD > simulation. I were basically referred to *>* the gromacs manual ( *>* > http://www.gromacs.org/Documentation/How-tos/Normal_Mode_Analysis > <http://www.gromacs.org/Documentation/How-tos/Normal_Mode_Analysis>) and > some *>* threads in the mailing list. *> >* First, Gromacs (and in some > literature) suggest to have L-BFGS energy *>* minimization (sometimes > after steep and cg calculations) of the structure *>* before NMA and > L-BFGS generally requires shifted/switched interactions. *>* However, for > vdW the option should have been changed by introducing the *>* > vdw-modifier and l-bfgs seems to still recognize it as cutoff. *> >* > Second, for electrostatics, the option has been improved by *>* > Reaction-field-zero but is it still suitable for NMA calculation? *> >* > Third, I have achieved a maximum force of 9.78750e-01 using cg and still * > >* got the warning: "The force is probably not small enough to ensure > that you *>* are at a minimum." Usually how stable do we have to achieve > in order to *>* avoid the warning? *> >* Thanks in advance for any > experience shared. *You need to compile gromacs in double precision to > get to low forces. > >* Kevin *>* OSU Pharmacy *> -- > David van der Spoel, Ph.D., Professor of Biology > Head of Department, Cell & Molecular Biology, Uppsala University. > Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205. > http://www.icm.uu.se -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Domain decomposition for parallel simulations
Dear Users, I have encountered the problem of "There is no domain decomposition for n nodes that is compatible with the given box and a minimum cell size of x nm" and by reading through the gromacs website and some threads I understand that the problem might be caused by breaking the system into too small boxes by too many ranks. However, I have no idea how to get the correct estimation of suitable paralleling parameters. Hope someone could share his experience. Here are information stated in the log file: *Initializing Domain Decomposition on 4000 ranks* *Dynamic load balancing: on* *Will sort the charge groups at every domain (re)decomposition* *Initial maximum inter charge-group distances:* *two-body bonded interactions: 0.665 nm, Dis. Rest., atoms 23558 23590* * multi-body bonded interactions: 0.425 nm, Proper Dih., atoms 12991 12999* *Minimum cell size due to bonded interactions: 0.468 nm* *Maximum distance for 5 constraints, at 120 deg. angles, all-trans: 0.819 nm* *Estimated maximum distance required for P-LINCS: 0.819 nm* *This distance will limit the DD cell size, you can override this with -rcon* *Guess for relative PME load: 0.11* *Will use 3500 particle-particle and 500 PME only ranks* *This is a guess, check the performance at the end of the log file* *Using 500 separate PME ranks, as guessed by mdrun* *Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25* *Optimizing the DD grid for 3500 cells with a minimum initial size of 1.024 nm* *The maximum allowed number of cells is: X 17 Y 17 Z 17* And I got this afterwards: *Fatal error:* *There is no domain decomposition for 3500 ranks that is compatible with the given box and a minimum cell size of 1.02425 nm* Here are some questions: 1. the maximum allowed number of cells is 17x17x17 which is 4913 and seems to be larger than the requested 3500 particle-particle ranks, so the minimum cell size is not causing the problem? 2. Where does this 1.024 nm comes from? We can see the inter charge-group distances are listed as 0.665 and 0.425 nm 3. The distance restraint between atoms 23558 23590 was set explicitly (or added manually) in the topology file and should be around 0.32 nm by using [intermolecular_interactions]. How could I know my manual setting is working or not? As it has shown a different value. Thanks in advance, Kevin OSU -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Precision of trr and checkpoint files
Dear Users, I am wondering how precision of the files affects the continuation of simulations. I have compiled both single and double precision of Gromacs 5.0.6 but I used to work with the double precision one. Usually I continue a simulation from a equilibration run (probably with Berendsen schemes) to a production run (with parrinello-rahman). So I have to use grommp instead of a tpbconv or -cpi option. I used to specify -c, -p and that's all. I have checked my gro file and it contains corrdinates, velocities and box size so it should be fine. However when I learnt from the command manual and used -t previous.trr to feed the corrdinates, velocities and box size, grommp just asked me for a conf.gro. Shouldn't a trr file contains all these things? Okay I still can go on with -c, -p plus -t (Though it is kind of duplication for me), but this time the precision of the files is not understandable. When I do not use -t, everything seems perfectly fine. When I use -t previous.cpt, it gives Reading Coordinates, Velocities and Box size from old trajectory Will read whole trajectory Precision mismatch for state entry box, code precision is double, file precision is float Precision mismatch for state entry box-rel, code precision is double, file precision is float Precision mismatch for state entry pres_prev, code precision is double, file precision is float Precision mismatch for state entry x, code precision is double, file precision is float Precision mismatch for state entry v, code precision is double, file precision is float Last frame -1 time 1000.000 Using frame at t = 1000 ps Starting time for run is 0 ps When I use -t previous.trr, it gives Reading Coordinates, Velocities and Box size from old trajectory Will read whole trajectory trn version: GMX_trn_file (single precision) Last frame 50 time 1000.000 Using frame at t = 1000 ps Starting time for run is 0 ps I think both screen-out means neither the trajectory nor the cpt file is in double precision, given that I have been used a double-precision Gromacs nearly all the time. I have gone through the online manual and some mailing threads but little has straightly answered my question. Also although it is stated on the manual ( http://www.gromacs.org/Documentation/Terminology/Precision?highlight=precision) and FAQ that single (mixed) precision of Gromacs is enough for production, can anyone elaborate a little bit more on the accuracy of this mixed precision to (usually very large) biological system simulations. As guys in our group used to run with (at least namely) double precision codes but the GPU speed-up was too attractive and Gromacs only provides single precision compiles for that. Thanks in advance, Kevin City University of Hong Kong -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Energy minimisation goes to several values
Thanks Justin again for the reply. Your example was very clear and I would rather drop such minimization protocol. Honestly such protocol comes from the conventional step-wise relaxation (or pre-equilibrium) practice for protein MD simulations. I believe you must have heard people constraining heavy atoms (carbon alpha atoms or sometimes even the backbone) during heating, box relaxation and then releasing them by decreasing the spring constant step-by-step before they started a production on the protein. So how do you comment on this stepwise protocol? If it is applicable, should it be turned on during energy minimization or just started from heating? Thanks in advance, Kevin On 27 Jun, 2015, at 07:54, gromacs.org_gmx-users-requ...@maillist.sys.kth.se gromacs.org_gmx-users-requ...@maillist.sys.kth.se wrote: From: Justin Lemkul jalem...@vt.edu mailto:jalem...@vt.edu Subject: Re: [gmx-users] Energy minimisation goes to several values Date: 27 June, 2015 07:54:13 HKT To: Discussion list for GROMACS users gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org Reply-To: gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org On 6/26/15 11:54 AM, Kevin C Chan wrote: Thanks Justin for the reply. Honestly I thought minimizing a large and complex (usually biomolecular) system part by part instead of as a whole will effectively shorten the computational time cost while causing no effect on the final structure. When you say “impedes”, do you mean it causes a longer calculation time in total or it will give a bad final structure? Well, you're running three minimizations instead of one, and you're achieving an unstable result. I'd say the three-step approach is not worth doing :) Consider something really simple - a polar, surface residue on a protein surrounded by water. Let's say you freeze the protein and let the water relax. The local waters respond to the fixed geometry of the side chain, which is (maybe) from a crystal and therefore perhaps not the correct conformation in solution. So the waters reorganize a bit. Then you let the protein relax but the waters are fixed. The side chain responds to a fixed clathrate of water that have been minimized around the wrong side chain conformation. What have you achieved? Nothing. Sure, you then minimize the whole system, but your starting point is potentially less plausible than it was to start with! At minimum, it's just a waste of time. Occasional use of restraints can be beneficial in some cases. Any time you talk about absolutely fixing large groups of atoms (like immobilizing water), I think it's really a waste of time. If a single, unrestrained minimization still leads to an unstable system, then it's not your minimization protocol that's to blame, rather an unresolvable starting structure or a bad topology. -Justin Thanks in advance, Kevin On 26 Jun, 2015, at 22:21, gromacs.org_gmx-users-requ...@maillist.sys.kth.se mailto:gromacs.org_gmx-users-requ...@maillist.sys.kth.se mailto:gromacs.org_gmx-users-requ...@maillist.sys.kth.se mailto:gromacs.org_gmx-users-requ...@maillist.sys.kth.se gromacs.org_gmx-users-requ...@maillist.sys.kth.se mailto:gromacs.org_gmx-users-requ...@maillist.sys.kth.se mailto:gromacs.org_gmx-users-requ...@maillist.sys.kth.se mailto:gromacs.org_gmx-users-requ...@maillist.sys.kth.se wrote: *From:*Justin Lemkul jalem...@vt.edu mailto:jalem...@vt.edu mailto:jalem...@vt.edu mailto:jalem...@vt.edu *Subject:**Re: [gmx-users] Energy minimisation goes to several values* *Date:*26 June, 2015 21:27:32 HKT *To:*gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org *Reply-To:*gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org On 6/25/15 9:45 PM, Kevin C Chan wrote: Dear Users, I am energy minimising a quite large solvated system containing protein and lipids (~800,000 atoms). I used to fix components of the system in order to speed-up energy minimisation and sometimes it is easier to debug such processes. Here is my protocol: 1. fix all except water and so to minimise water 2. fix water and then minimise all the rest atoms 3. fix nothin and then minimise the whole system While monitoring the energy of the system thought minimisations, it goes fine for step 1 and 2 and converged after just few hundred steps. However it goes back to several higher values of energy (bouncing between the values) and they started to increase very slowly for step 3. This makes no sense to me and did anyone have a similar experience? There are two unusual points: 1. The system energy drops suddenly instead of decreased gradually during step2 and then stays at a constant value. 2. If I use the resulting structure from step3 to proceed a, say, heating process, it simply blows up. To be clear, my system was solvated
Re: [gmx-users] Energy minimisation goes to several values
Thanks Justin for the reply. Honestly I thought minimizing a large and complex (usually biomolecular) system part by part instead of as a whole will effectively shorten the computational time cost while causing no effect on the final structure. When you say “impedes”, do you mean it causes a longer calculation time in total or it will give a bad final structure? Thanks in advance, Kevin On 26 Jun, 2015, at 22:21, gromacs.org_gmx-users-requ...@maillist.sys.kth.se gromacs.org_gmx-users-requ...@maillist.sys.kth.se wrote: From: Justin Lemkul jalem...@vt.edu mailto:jalem...@vt.edu Subject: Re: [gmx-users] Energy minimisation goes to several values Date: 26 June, 2015 21:27:32 HKT To: gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org Reply-To: gmx-us...@gromacs.org mailto:gmx-us...@gromacs.org On 6/25/15 9:45 PM, Kevin C Chan wrote: Dear Users, I am energy minimising a quite large solvated system containing protein and lipids (~800,000 atoms). I used to fix components of the system in order to speed-up energy minimisation and sometimes it is easier to debug such processes. Here is my protocol: 1. fix all except water and so to minimise water 2. fix water and then minimise all the rest atoms 3. fix nothin and then minimise the whole system While monitoring the energy of the system thought minimisations, it goes fine for step 1 and 2 and converged after just few hundred steps. However it goes back to several higher values of energy (bouncing between the values) and they started to increase very slowly for step 3. This makes no sense to me and did anyone have a similar experience? There are two unusual points: 1. The system energy drops suddenly instead of decreased gradually during step2 and then stays at a constant value. 2. If I use the resulting structure from step3 to proceed a, say, heating process, it simply blows up. To be clear, my system was solvated and auto-ionized using VMD tools and some water inside the membrane has been directly deleted. Backbone of the protein and phosphorus atoms of the membrane are under a position constraint during all the minimisations. I was choosing conjugate gradient for minimization. Does a normal minimization (just one overall minimization with nothing fixed) yield a stable starting point? Fixing atoms (using freezegrps?) often actually impedes minimization. -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu mailto:jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Energy minimisation goes to several values
Dear Users, I am energy minimising a quite large solvated system containing protein and lipids (~800,000 atoms). I used to fix components of the system in order to speed-up energy minimisation and sometimes it is easier to debug such processes. Here is my protocol: 1. fix all except water and so to minimise water 2. fix water and then minimise all the rest atoms 3. fix nothin and then minimise the whole system While monitoring the energy of the system thought minimisations, it goes fine for step 1 and 2 and converged after just few hundred steps. However it goes back to several higher values of energy (bouncing between the values) and they started to increase very slowly for step 3. This makes no sense to me and did anyone have a similar experience? There are two unusual points: 1. The system energy drops suddenly instead of decreased gradually during step2 and then stays at a constant value. 2. If I use the resulting structure from step3 to proceed a, say, heating process, it simply blows up. To be clear, my system was solvated and auto-ionized using VMD tools and some water inside the membrane has been directly deleted. Backbone of the protein and phosphorus atoms of the membrane are under a position constraint during all the minimisations. I was choosing conjugate gradient for minimization. Thanks in advance, Kevin City University of Hong Kong -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] position restraints
Sorry for digging the old thread, but I really have question concerning Justin’s reply. On 17 May, 2015, at 05:05, Justin Lemkul jalem...@vt.edu wrote: On 5/16/15 8:13 AM, Ming Tang wrote: Dear all, I want to minimize my system to relax the solvent with the positions of heavy atoms within the protein restrained. I searched the archives and got the idea that I need to put something like this in topology: #ifdef POSRES_WATER ; Position restraint for each water oxygen [ position_restraints ] ; i funct fcxfcyfcz 11 1000 1000 1000 #endif But I still have no idea how to get the heavy atoms restrained. Could anybody tell me how to do this? pdb2gmx does this automatically for you. It writes posre.itp, which specifies heavy atom restraints. It's controlled by a similar #ifdef block. I want to use gromos54a7 ff to simulate my collagen. As Justin suggested, I need to study the original paper Definition and testing of the GROMOS force-field versions 54A7 and 54B7 to write my .mdp file. In the paper it says, 1. Each system was energy minimised to relax the solvent with the positions of the heavy atoms within the protein and RNA restrained. 2. All solute atom positions were restrained to their positions in the initial structure through a harmonic potential energy term with force constant of 2.5 9 104 kJ mol-1 nm-2. Could anybody tell me how to restrain the heavy atoms of protein, and how to get the solute atom positions restrained with a harmonic potential energy term with a certain constant? This is what posre.itp does. Also note that you don't necessarily have to follow that protocol exactly; my general advice for proper use of force fields is about proper nonbonded setup. Preparation of a system often depends on what the system is and if there are any special considerations. For a simple protein system, it's pretty hard to break things. I don't see much benefit in restraining atoms during minimization, though. I believe a significant number of the users learn MD protocols from reading papers and often it is difficult to trace back the reasons behind certain steps. Restraining atoms during minimization or even heating-up, pre-equilibrium is a very good example. While Justin mentioned that the step is kind of “useless”, I read quite a amount of similar approaches in papers that people “relax” their system step by step before a production. I will greatly appreciate if anyone can provide some original papers on discussing rigorous MD protocols (especially for nowadays complex biomolecular systems). Thanks in advance, -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. Kevin PhD Candidate Department of Physics and Material Science City University of Hong Kong uk...@gmx.hk -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Changing selections in Umbrella Sampling
Dear Users, I am doing umbrella sampling of pulling a protein away from a membrane. Initially I chose reaction coordinate to be the distance between com of the lower layer of the lipids and of the protein. After weeks of simulations, I have concluded a better selection of com - as the protein is positioned away from the center, only the membrane underneath should be used to measure the distance. Should I start from the very beginning, or I can actually stick to the last frame of the previous umbrella sampling and save times for equilibration? Looking forward to any experience sharing. Thanks in advance, Kevin City University of Hong Kong -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] No default Ryckaert-Bell. types for proline
Thanks Justin, I'm sorry that I'm not clear with what do you mean by so you can likely just use the same parameters in ffbonded.itp for this dihedral in the zwitterionic form. Also as I have mentioned before, when I look into the line 161, it refers to a dihedral term concerning (CD N CA C) which is not able to be found in ffbonded.itp. However, other dihedral terms such as (CD N CA CB) or (CD N CA HA) which are also not found in itp do not report such error and I was successful to grompp after commenting out the line 161. So why would this happen? Thanks in advance, Kevin On 2 Jun, 2015, at 20:05, Justin Lemkul jalem...@vt.edu javascript:_e(%7B%7D,'cvml','jalem...@vt.edu'); wrote: On 6/1/15 12:25 PM, Kevin C Chan wrote: pdb2gmx command: gmx_d pdb2gmx -f proline.pdb -ter 15: (OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)) 1: TIP4P TIP 4-point, recommended 1: PRO-ZWITTERION_NH2+ (only use with zwitterions containing exactly one residue) 1: PRO-ZWITTERION_COO- (only use with zwitterions containing exactly one residue) (I have also tried 4: ZWITTERION_NH3+ (only use with zwitterions containing exactly one residue) but the result was the same) Screen-out: http://pastebin.com/CuMa7aZ7# topology: http://pastebin.com/FSbcFGYR It appears that there are simply no parameters in the OPLS-AA force field for this entity. Whether or not this is because OPLS-AA doesn't officially have parameters for it or they are just missing from the GROMACS implementation is something you'd have to dig into, but largely it seems that all dihedral parameters for the N- and C-terminal forms of PRO are all the same, so you can likely just use the same parameters in ffbonded.itp for this dihedral in the zwitterionic form. -Justin Thanks in advance, Kevin On Tue, Jun 2, 2015 at 12:11 AM, Justin Lemkul jalem...@vt.edu javascript:_e(%7B%7D,'cvml','jalem...@vt.edu'); wrote: On 6/1/15 12:09 PM, Kevin C Chan wrote: Thanks Justin, I did not but I tried just now. The aforementioned dihedral terms remain along with the errors. Please provide your full pdb2gmx command, all screen output, and full topology. Use file sharing services, pastebin, etc. as needed. -Justin Thanks in advance, Kevin On Mon, Jun 1, 2015 at 11:13 PM, Justin Lemkul jalem...@vt.edu javascript:_e(%7B%7D,'cvml','jalem...@vt.edu'); wrote: On 6/1/15 10:57 AM, Kevin C Chan wrote: Dear Users, I am testing with OPLSSaa force-field of proline. First I arbitrarily extract a single PRO from a pdb and fed into pdb2gmx. The resulted gro and top are however not suitable for further grompp (even minimization). The error messages are as follows: ERROR 1 [file topol.top, line 161]: No default Ryckaert-Bell. types When I look into the line 161, it refers to a dihedral term concerning (CD N CA C) which is not able to be found in ffbonded.itp. However, other dihedral terms such as (CD N CA CB) or (CD N CA HA) which are also not found in itp do not report such error and I was successful to grompp after commenting out the line 161. Why would this happen? Did you choose the proper PRO-ZWITTERION-* termini when running pdb2gmx? -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu javascript:_e(%7B%7D,'cvml','jalem...@outerbanks.umaryland.edu'); | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org javascript:_e(%7B%7D,'cvml','gmx-users-requ...@gromacs.org');. -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu javascript:_e(%7B%7D,'cvml','jalem...@outerbanks.umaryland.edu'); | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx
Re: [gmx-users] No default Ryckaert-Bell. types for proline
Thanks Justin, I have long been confused by atom names and types in Gromacs topology. I have just went through the manual again, and it confirms me that what we read directly in ffbonded.itp under GMXHOME/share/gromacs/top/oplsaa.ff are atomtypes. In other words, what you have listed in the last email C CT_2 N CT_33 -5.72371 -18.33847 -5.23419 29.29636 0.0 0.0 ; Pro (fit to AM1) CD-N-CA-C, JT-R 2/10/97 C CT_2 NT CT_33 -5.72371 -18.33847 -5.23419 29.29636 0.0 0.0 ; Pro (fit to AM1) CD-N-CA-C, JT-R 2/10/97 C_3CT_2 N CT_33 -5.72371 -18.33847 -5.23419 29.29636 0.0 0.0 ; Pro COO- terminus. C_3CT_2 NT CT_33 -5.72371 -18.33847 -5.23419 29.29636 0.0 0.0 ; Pro COO- terminus. which was extracted from ffbonded.itp is composed of atomtypes. Meanwhile what I have stated in the brackets - (CD N CA C), (CD N CA CB) or (CD N CA HA) are by myself read from the topology file after pdb2gmx. To be exact, I read [ dihedrals ] ; aiajakal functc0c1c2 c3c4c5 12 1 4 5 3 12 1 4 6 3 12 1 415 3 from the topology and refer the numbers (e.g. 12 1 4 15) to the above directive [ atoms ] to find out that it is (CD N CA C). We can also read from directive [ atoms ] that under the type column, it should be (opls_296 opls_309 opls_285 opls_271). So which is the atomtypes, (CD N CA C) or (opls_296 opls_309 opls_285 opls_271)? It seems (opls_296 opls_309 opls_285 opls_271) are the atomtypes, but why the itp states (C CT_2 N CT_3)? And how can I refer my topology to this line of parameter when it deals (CD N CA C) or (opls_296 opls_309 opls_285 opls_271)? Regards, Kevin On Wed, Jun 3, 2015 at 10:11 AM, Justin Lemkul jalem...@vt.edu wrote: On 6/2/15 8:51 PM, Kevin C Chan wrote: Thanks Justin, I'm sorry that I'm not clear with what do you mean by so you can likely just use the same parameters in ffbonded.itp for this dihedral in the zwitterionic form. The Pro dihedral term for this torsion is the same, irrespective of whether it's N- or C-terminal. From ffbonded.itp: C CT_2 N CT_33 -5.72371 -18.33847 -5.23419 29.29636 0.0 0.0 ; Pro (fit to AM1) CD-N-CA-C, JT-R 2/10/97 C CT_2 NT CT_33 -5.72371 -18.33847 -5.23419 29.29636 0.0 0.0 ; Pro (fit to AM1) CD-N-CA-C, JT-R 2/10/97 C_3CT_2 N CT_33 -5.72371 -18.33847 -5.23419 29.29636 0.0 0.0 ; Pro COO- terminus. C_3CT_2 NT CT_33 -5.72371 -18.33847 -5.23419 29.29636 0.0 0.0 ; Pro COO- terminus. All the same, so the parameters for the zwitterionic species can likely be copied to your missing interaction. Also as I have mentioned before, when I look into the line 161, it refers to a dihedral term concerning (CD N CA C) which is not able to be found in ffbonded.itp. However, other dihedral terms such as (CD N CA CB) or (CD N CA HA) which are also not found in itp do not report such error and I was Atom names and atom types are entirely different concepts. grompp doesn't care what you call your atoms (atom name), but it does care what chemical characteristics they have (atom type). successful to grompp after commenting out the line 161. So why would this happen? Don't ever comment out lines because of issues like this. By doing so, you're saying that torsion contributes nothing to the energy function, which is wrong. I suggest copying the above parameters to the correct sequence of atom types. Stuff like this happens all the time. Force fields are massive entities, and are added to incrementally as need arises. A zwitterionic proline is probably not something that people commonly use, and its the most complex case due to the chemical nature of it. It's entirely possible that the .tdb entries were simply added by analogy to existing groups as a convenience and never fully tested. -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive
Re: [gmx-users] No default Ryckaert-Bell. types for proline
Thanks Justin, I did not but I tried just now. The aforementioned dihedral terms remain along with the errors. Thanks in advance, Kevin On Mon, Jun 1, 2015 at 11:13 PM, Justin Lemkul jalem...@vt.edu wrote: On 6/1/15 10:57 AM, Kevin C Chan wrote: Dear Users, I am testing with OPLSSaa force-field of proline. First I arbitrarily extract a single PRO from a pdb and fed into pdb2gmx. The resulted gro and top are however not suitable for further grompp (even minimization). The error messages are as follows: ERROR 1 [file topol.top, line 161]: No default Ryckaert-Bell. types When I look into the line 161, it refers to a dihedral term concerning (CD N CA C) which is not able to be found in ffbonded.itp. However, other dihedral terms such as (CD N CA CB) or (CD N CA HA) which are also not found in itp do not report such error and I was successful to grompp after commenting out the line 161. Why would this happen? Did you choose the proper PRO-ZWITTERION-* termini when running pdb2gmx? -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] No default Ryckaert-Bell. types for proline
pdb2gmx command: gmx_d pdb2gmx -f proline.pdb -ter 15: (OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)) 1: TIP4P TIP 4-point, recommended 1: PRO-ZWITTERION_NH2+ (only use with zwitterions containing exactly one residue) 1: PRO-ZWITTERION_COO- (only use with zwitterions containing exactly one residue) (I have also tried 4: ZWITTERION_NH3+ (only use with zwitterions containing exactly one residue) but the result was the same) Screen-out: http://pastebin.com/CuMa7aZ7# topology: http://pastebin.com/FSbcFGYR Thanks in advance, Kevin On Tue, Jun 2, 2015 at 12:11 AM, Justin Lemkul jalem...@vt.edu wrote: On 6/1/15 12:09 PM, Kevin C Chan wrote: Thanks Justin, I did not but I tried just now. The aforementioned dihedral terms remain along with the errors. Please provide your full pdb2gmx command, all screen output, and full topology. Use file sharing services, pastebin, etc. as needed. -Justin Thanks in advance, Kevin On Mon, Jun 1, 2015 at 11:13 PM, Justin Lemkul jalem...@vt.edu wrote: On 6/1/15 10:57 AM, Kevin C Chan wrote: Dear Users, I am testing with OPLSSaa force-field of proline. First I arbitrarily extract a single PRO from a pdb and fed into pdb2gmx. The resulted gro and top are however not suitable for further grompp (even minimization). The error messages are as follows: ERROR 1 [file topol.top, line 161]: No default Ryckaert-Bell. types When I look into the line 161, it refers to a dihedral term concerning (CD N CA C) which is not able to be found in ffbonded.itp. However, other dihedral terms such as (CD N CA CB) or (CD N CA HA) which are also not found in itp do not report such error and I was successful to grompp after commenting out the line 161. Why would this happen? Did you choose the proper PRO-ZWITTERION-* termini when running pdb2gmx? -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] No default Ryckaert-Bell. types for proline
Dear Users, I am testing with OPLSSaa force-field of proline. First I arbitrarily extract a single PRO from a pdb and fed into pdb2gmx. The resulted gro and top are however not suitable for further grompp (even minimization). The error messages are as follows: ERROR 1 [file topol.top, line 161]: No default Ryckaert-Bell. types When I look into the line 161, it refers to a dihedral term concerning (CD N CA C) which is not able to be found in ffbonded.itp. However, other dihedral terms such as (CD N CA CB) or (CD N CA HA) which are also not found in itp do not report such error and I was successful to grompp after commenting out the line 161. Why would this happen? Thanks in advance, Kevin City University of Hong Kong -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Is it reasonable to get a bad Ramachandran plot after MD?
Dear Justin, Thanks for the reply. I was using CHARMM27 force-field. In fact, it was a MDFF simulations in NAMD raising my interest as the resulted structure had 5% outliers and was doubted by experimentalists. I also checked other conventional MD results and the % of outliers vary among structures from 1% to 5% (all using the same CHARMM force-field). I do not know specifically which residues lie in those regions, can anyone suggest a more comprehensive analysis tool for Ramachandran plot? Thanks in advance, Kevin On Fri, May 15, 2015 at 12:41 AM, Justin Lemkul jalem...@vt.edu wrote: On 5/14/15 8:51 AM, Kevin C Chan wrote: Dear Users, As the title states, I am wondering the validity behind the outliers appear in a Ramachandran plot of structure after MD. I only have little data on this, but for a 10,000-atom protein, portion of outlier can increase from 0.7% to 5.0%. Is this explainable? Or it simply means my force-field does not care about the dihedral angles but still accurate enough to predict protein dynamics? Which residues lie in those regions? What force field are you using? Is that force field known to produce only-allowable Ramachandran regions. Note that glycine can do just about anything it likes... -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] (no subject)
Dear Users, I'm currently doing analysis of results of umbrella sampling. To guarantee the reliability of the PMF profile, I plot the histograms to check the overlapping of them. While we have to set a bin size for the analysis, I found the bin size affect the overlapping of the histograms. For instance, a larger bin size (e.g. 1600 for a PMF profile over 30A distance) will cause a visually poorer overlapping of the histograms than a smaller bin size (e.g. 200). 1) Is this normal? Can anyone share any experience on this? 2) What is the best bin size for an umbrella sampling analysis? Any answer is highly appreciated. Thanks in advance, Kevin City University of Hong Kong -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Bin size for Umbrella Sampling Analysis
Dear Users, I'm currently doing analysis of results of umbrella sampling. To guarantee the reliability of the PMF profile, I plot the histograms to check the overlapping of them. While we have to set a bin size for the analysis, I found the bin size affect the overlapping of the histograms. For instance, a larger bin size (e.g. 1600 for a PMF profile over 30A distance) will cause a visually poorer overlapping of the histograms than a smaller bin size (e.g. 200). 1) Is this normal? Can anyone share any experience on this? 2) What is the best bin size for an umbrella sampling analysis? -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Bin size for an Umbrella Sampling analysis
Dear Users, I am doing analysis for umbrella sampling and find the chosen bin size affecting the overlapping of histograms. While we should always guarantee a “well” overlapped histogram, choosing of bin size seems to have certain criteria that I did not expect before. Looking forward to any experience sharing. Regards, Kevin PhD Candidate Department of Physics and Material Science City University of Hong Kong uk...@gmx.hk -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Langevin temperature coupling in Gromacs
Dear Users, I am thinking whether a Langevin thermostat can be achieved by turning on integrator=sd in Gromacs and setting bd-fric. All I want is to make an equivalent mdp file (thermostat part) compared to an NAMD parameter in which the Langevin thermostat is turned on and gamma equals to 5. Thanks in advance, Kevin -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] gromacs.org_gmx-users Digest, Vol 128, Issue 91
Thanks for the reply. In principle that is what I plan to do - generate a basic topology (with homogeneous charge information) and then manually (or a script) change the charge of each SI atom. Honestly I only have experience of writing a simple topology for a MARTINI system. So what I expect was first to include the gromos forcefield so I have the [defaults] [atomtypes] [bondtypes] defined and then manually write the [moleculetypes] and [atoms] in which I can define the charges (I assume it will overwrite the charge value defined in [atomtypes], right?) then the tens or hundreds of rows of entries were treated as one single molecule and finally the [system] and [molecules]. However I am now stuck at the step that it seems no forcefield provided by Gromacs gives bonded parameters for Si-Si so g_x2top could not work out. Or maybe I should give up g_x2top and just manually write the [bondtypes] according some parameters provided by the literature and combined with nonbonded parameters provided by gromos forcefield, for instance, but this then arises the consistency problem Oh god :( Thanks in advance and looking forward to any comment or experience sharing. Regards, Kevin uk...@gmx.hk From: Justin Lemkul jalem...@vt.edu To: gmx-us...@gromacs.org Date: 20 December, 2014 06:11:14 HKT Reply-To: gmx-us...@gromacs.org Subject: Re: [gmx-users] Preparing a silicon bulk On 12/19/14 12:43 PM, Kevin C Chan wrote: I am kind of confused, how could I actually reply and include my reply into a thread in the mailing list? As I am receiving digest email, I should not just press the reply bottom, or should I? Replies work like any other email. Click reply, change the subject, chop out the text that doesn't pertain to your issue, and send. If digests are acting weird, switch to individual messages. The list can be somewhat high-traffic, but certainly the last few weeks haven't been all that busy. Back to the question, thank you for providing the hint on unit of Gromacs, g_x2top could identify the bonds now. However g_x2top is still not working as: Opening force field file /home/kevin/opt/gmx_fftw3_double/share/gromacs/top/gromos54a7.ff/atomname2type.n2t There are 1 name to type translations in file gromos54a7.ff Generating bonds from distances... atom 8 Can not find forcefield for atom SI1-1 with 4 bonds Can not find forcefield for atom SI2-2 with 4 bonds Can not find forcefield for atom SI3-3 with 4 bonds Can not find forcefield for atom SI4-4 with 4 bonds Can not find forcefield for atom SI5-5 with 4 bonds Can not find forcefield for atom SI6-6 with 4 bonds Can not find forcefield for atom SI7-7 with 4 bonds Can not find forcefield for atom SI8-8 with 4 bonds And my atonname2type.n2t contains: SI SI 1 28.08 3SI 0.2352 SI 0.2352 SI 0.2352 We could see that it identifies 4 bonds for each atom (but why not 3?), however no ff could be found. I look into gromos54a7 (I thought it contains Si forcefields) and I found no Si-Si bonded parameters in ffbonded.itp but a comprehensive set of non-bonded parameters of Si in ffnonbonded.itp. I am curious why it is like this. Did people only do single Si simulations? Si should form 4 bonds; that's what satisfies its valence. What g_x2top is telling you is: I found 4 atoms within the bond distance, but you tell me in the .n2t file that there should only be 3. The Si parameters in the force field were probably for silica, so no Si-Si bonds. Ultimately is it better (or easier) to manually construct the topology file? In theory, you have a regular ordering of atoms on a lattice, so g_x2top is the perfect tool, but the desire to assign different charges makes that difficult. I would suggest building the topology with g_x2top and modifying as needed; that will require some scripting (or could be done by hand for a reasonably small system), but realize that changes in charges can affect the validity of bonded parameters, since everything in a force field is interconnected. So even if you find suitable bonded parameters for the interactions you need, if you change aspects of the nonbonded parameters, the bondeds can be compromised as a result... -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit
Re: [gmx-users] Preparing a silicon bulk
Thanks for the reply. I am quite new to the field of Materials so I do need someone else to find a “good” paper on parameters for a silicon surface. While I am waiting to for that, I am actually very curious that why gromos does not have such parameters and nearly all popular forcefields (for biomolecule simulations) contains silica parameters. It is strange that we even have Si-Si nonbonded parameters in gromos but no bonded ones. In practice, it is also hard to say how my charges vary across the surface (they are just 0.X in difference). Basically I am thinking if I could generate a topology for a simple silicon bulk then I would just rewrite the [atoms] in [moleculetypes] to assign the charges. Like the follows: [ defaults ] ; nbfunc comb-rulegen-pairs fudgeLJ fudgeQQ 1 1no1.0 1.0 [ atomtypes ] ;name at.num mass charge ptype c6 c12 SI 1 28.08 0.000 A 1.000 1.000 [ bondtypes ] ;i j func b0 kb SI SI 1 1.000 1.000 #include gromos54a7.ff/forcefield.itp [ moleculetype ] ; molname nrexcl SIB 4 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 SI 1 SIB SI 1 -0.10 2 SI 1 SIB SI 2 0.60 3 SI 1 SIB SI 3 -0.30 4 SI 1 SIB SI 4 -0.20 [ system ] SIB in Water [ molecules ] SIB 1 In this way what I have to do is to plug correct numbers into [bondtypes]. Thanks in advance, Kevin uk...@gmx.hk On 20 Dec, 2014, at 22:13, Justin Lemkul jalem...@vt.edu wrote: On 12/20/14 5:29 AM, Kevin C Chan wrote: Thanks for the reply. In principle that is what I plan to do - generate a basic topology (with homogeneous charge information) and then manually (or a script) change the charge of each SI atom. Honestly I only have experience of writing a simple topology for a MARTINI system. So what I expect was first to include the gromos forcefield so I have the [defaults] [atomtypes] [bondtypes] defined and then manually write the [moleculetypes] and [atoms] in which I can define the charges (I assume it will overwrite the charge value defined in [atomtypes], right?) then the tens or hundreds of rows of entries were treated as one single molecule and finally the [system] and [molecules]. The charge values in ffnonbonded.itp are never used for anything; they're an artifact of old intentions. However I am now stuck at the step that it seems no forcefield provided by Gromacs gives bonded parameters for Si-Si so g_x2top could not work out. Or maybe I should give up g_x2top and just manually write the [bondtypes] according some parameters provided by the literature and combined with nonbonded parameters provided by gromos forcefield, for instance, but this then arises the consistency problem Oh god :( Find a suitable force field in the literature, something that you can use to reproduce their findings; implement it (something like an Si or silica surface should be very easy, given the few parameters that are needed) and test that you can make a sane topology. Then do your hacks and see how things go. I don't know what kind of validation you'd need to do, and you haven't said by how much, or why, the charges will vary across the surface, so that's also hard to say. Welcome to the delightful world of force field parametrization :) -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Preparing a silicon bulk
I am kind of confused, how could I actually reply and include my reply into a thread in the mailing list? As I am receiving digest email, I should not just press the reply bottom, or should I? Back to the question, thank you for providing the hint on unit of Gromacs, g_x2top could identify the bonds now. However g_x2top is still not working as: Opening force field file /home/kevin/opt/gmx_fftw3_double/share/gromacs/top/gromos54a7.ff/atomname2type.n2t There are 1 name to type translations in file gromos54a7.ff Generating bonds from distances... atom 8 Can not find forcefield for atom SI1-1 with 4 bonds Can not find forcefield for atom SI2-2 with 4 bonds Can not find forcefield for atom SI3-3 with 4 bonds Can not find forcefield for atom SI4-4 with 4 bonds Can not find forcefield for atom SI5-5 with 4 bonds Can not find forcefield for atom SI6-6 with 4 bonds Can not find forcefield for atom SI7-7 with 4 bonds Can not find forcefield for atom SI8-8 with 4 bonds And my atonname2type.n2t contains: SI SI 1 28.08 3SI 0.2352 SI 0.2352 SI 0.2352 We could see that it identifies 4 bonds for each atom (but why not 3?), however no ff could be found. I look into gromos54a7 (I thought it contains Si forcefields) and I found no Si-Si bonded parameters in ffbonded.itp but a comprehensive set of non-bonded parameters of Si in ffnonbonded.itp. I am curious why it is like this. Did people only do single Si simulations? Ultimately is it better (or easier) to manually construct the topology file? Regards, On 12/18/14 11:46 PM, Kevin C Chan wrote: Due to some reasons I could not receive digests from the list so I could only copy and paste from the archives and reply to you. I really hope this reply could go to the same thread of my question so it would not screw the list archive. Nope, this is the start of a new thread now. Make sure you sort out whatever delivery issues you're having using the admin page. Back to the question, first thank you so much for your reply Justin. I understand that pdb2gmx might not be useful, I just tried it. g_x2top seems better but I could not get it work even following the file format: it says it could find 1 name to type translations but it can not find forcefield for atoms like Opening force field file /home/kevin/opt/gmx_fftw3_double/share/gromacs/top/gromos54a7.ff/atomname2type.n2t There are 2 name to type translations in file gromos54a7.ff Generating bonds from distances... atom 8 Can not find forcefield for atom SI1-1 with 0 bonds Can not find forcefield for atom SI2-2 with 0 bonds Can not find forcefield for atom SI3-3 with 0 bonds Can not find forcefield for atom SI4-4 with 0 bonds Can not find forcefield for atom SI5-5 with 0 bonds Can not find forcefield for atom SI6-6 with 0 bonds Can not find forcefield for atom SI7-7 with 0 bonds Can not find forcefield for atom SI8-8 with 0 bonds I gave the n2t like SI SI 1 28.08 3SI 2.352 SI 2.352 SI 2.352 ;2.352 is measured from my input coordinates Gromacs uses nm for everything, not Angstrom. Given that bonds are detected based on a 10% distance tolerance, this is the reason for the above messages. and just trying with 8 SI atoms and in the future I definite would like to have more. While waiting for your comments on this g_x2top problem, I would also like to have suggestions on manually writing the top file. As I don't think we need those tedious dihedrals for a silicon bulk, I though the topology should look rather simple: #include gromos54a7.ff ; I have looked into gromos54a7 and found parameters for SI so this include is enough am I right? [ moleculetype ] ; name SomeName SI 8 [ atoms ] 1 SI 1 SI SI1 1 charge1 2 SI 1 SI SI1 1 charge2 ... 8 SI 1 SI SI1 1 charge8 ; Maybe bonds and angles are also needed? I am not sure You need to determine a suitable force field based on something that is published and validated to work. I don't know what the source of Si parameters is for Gromos96, but people use different force fields in published work. If you're wondering if bonds and angles are necessary (certainly the former, likely the latter, and maybe/maybe not on dihedrals, depending on the parametrization) then you need to back up and do some evaluation of parameters that are already available. The output is only as good as the input, and Gromacs tools won't work (too much) magic to save you from mistakes. -Justin #include spc.itp [ system ] SomeName [ molecules ] SI 8 SOL 100 Is this the case? Thanks in advance, Kevin Kevin uk...@gmx.hk -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send
[gmx-users] Preparing a silicon bulk
Due to some reasons I could not receive digests from the list so I could only copy and paste from the archives and reply to you. I really hope this reply could go to the same thread of my question so it would not screw the list archive. Back to the question, first thank you so much for your reply Justin. I understand that pdb2gmx might not be useful, I just tried it. g_x2top seems better but I could not get it work even following the file format: it says it could find 1 name to type translations but it can not find forcefield for atoms like Opening force field file /home/kevin/opt/gmx_fftw3_double/share/gromacs/top/gromos54a7.ff/atomname2type.n2t There are 2 name to type translations in file gromos54a7.ff Generating bonds from distances... atom 8 Can not find forcefield for atom SI1-1 with 0 bonds Can not find forcefield for atom SI2-2 with 0 bonds Can not find forcefield for atom SI3-3 with 0 bonds Can not find forcefield for atom SI4-4 with 0 bonds Can not find forcefield for atom SI5-5 with 0 bonds Can not find forcefield for atom SI6-6 with 0 bonds Can not find forcefield for atom SI7-7 with 0 bonds Can not find forcefield for atom SI8-8 with 0 bonds I gave the n2t like SI SI 1 28.08 3SI 2.352 SI 2.352 SI 2.352 ;2.352 is measured from my input coordinates and just trying with 8 SI atoms and in the future I definite would like to have more. While waiting for your comments on this g_x2top problem, I would also like to have suggestions on manually writing the top file. As I don't think we need those tedious dihedrals for a silicon bulk, I though the topology should look rather simple: #include gromos54a7.ff ; I have looked into gromos54a7 and found parameters for SI so this include is enough am I right? [ moleculetype ] ; name SomeName SI 8 [ atoms ] 1 SI 1 SI SI1 1 charge1 2 SI 1 SI SI1 1 charge2 ... 8 SI 1 SI SI1 1 charge8 ; Maybe bonds and angles are also needed? I am not sure #include spc.itp [ system ] SomeName [ molecules ] SI 8 SOL 100 Is this the case? Thanks in advance, Kevin On 12/18/14 12:38 PM, Kevin C Chan wrote: Dear Users, I wish to simulate a Silicon bulk using Gromacs. I have created the bulk using a plugin in VMD but I found difficulties when generating a top file from this pdb file. I have tried pdb2gmx and x2top but both does not work. Gromacs tells me that although we already have parameters for SI in maybe ffnonbonded.itp, we do not have SI in the residue topology database and also ffbonded.itp. How could it be and I do appreciate if anyone could share his experience of simulating a silicon bulk in Gromacs. pdb2gmx is not the answer, since it was designed for linear polymers with limited branching. The only option for materials like this is g_x2top. That requires creation of a suitable .n2t file, the format for which is described in http://www.gromacs.org/Documentation/File_Formats/.n2t_File - the contents should be fairly simple (one or two lines?) for silica. But... Another thing, I would like to assign different charge to different SI atom in my bulk. How can I achieve this? ...if you want to do this, there is no automated method. You could, in theory, assign different residue names to each Si unit and assign that charge, but then the parameters of O are likely to be contingent upon the type of Si to which it is attached. This eliminates g_x2top, which requires uniformity. And as said above, this is not a task suited for pdb2gmx. So either you need to (1) use g_x2top and hack or (2) write your own sort of code/script to generate the topology. -Justin -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Genbox question
Dear Users, I am trying to build some starting structures of lipid bilayers for MD using genbox to add waters. For some reasons I have already equilibrated my lipids in a water box previously so it no longer looks like a perfect rectangle, there are “flying tails at the edges. However I still want to build a typical bilayer-in-water model in which water only appears above and below the bilayer (but not asides) and it should be able to introduce a periodic boundary condition in the x-y plane (assuming the bilayer’s normal is along the z-axis) afterwards. The problem is as I could not accurately define the dimensions of my equilibrated lipids anymore, how could genbox be used to fill the box as what I expected? I will greatly appreciate it if anyone could share your experience on this. Cheers, Kevin -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
