Re: [gmx-users] Forcefield parameters
On 28/10/2010 3:03 AM, Sai Pooja wrote: Hi Mark, I am familiar with the section, however, I have a few doubts. When you say ... So modified protein-water VDW interactions can be introduced by defining all relevant protein atom-TIP3P oxygen [nonbond_params] terms. Do you mean that I can introduce [pairtypes] specifying interactions of all relevant protein atoms with tip3p atoms? No, [pairtypes] are not relevant. and I do not understand - It may be simpler to modify the [atomtypes] to generate the modified VDW from the combination rule Would this not alter all non-bonded interactions for the atomtypes defined? I think I haven't understood thsi properly... I was thinking of a flawed procedure. Ignore this part. Mark Pooja On Tue, Oct 26, 2010 at 8:21 PM, Mark Abraham mark.abra...@anu.edu.au mailto:mark.abra...@anu.edu.au wrote: I think manual 5.3.3 covers the relevant points. Let me know what is not clear. Mark - Original Message - From: Sai Pooja saipo...@gmail.com mailto:saipo...@gmail.com Date: Wednesday, October 27, 2010 9:41 Subject: Re: [gmx-users] Forcefield parameters To: Discussion list for GROMACS users gmx-users@gromacs.org mailto:gmx-users@gromacs.org On Tue, Oct 26, 2010 at 6:04 PM, Mark Abraham mark.abra...@anu.edu.au wrote: - Original Message - From: Sai Pooja saipo...@gmail.com Date: Wednesday, October 27, 2010 8:52 Subject: Re: [gmx-users] Forcefield parameters To: Discussion list for GROMACS users gmx-users@gromacs.org On Tue, Oct 26, 2010 at 4:04 PM, Justin A. Lemkul jalem...@vt.edu wrote: Sai Pooja wrote: On Tue, Oct 26, 2010 at 3:46 PM, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Sai Pooja wrote: Hi, I want to change the non-bonded parameters to modify the interaction between water molecules and protein molecules. I am using CHARMM forcefield with Tip3p water. The ffnonbonded.itp file of the forcefield has non-bonded parameters for tip3p water. Can I achieve the above by changing these parameters? That depends on your definition of modify, but yes, in a way, you can make changes here. 1) Modify - Multiply sigma and epsilon by a constant If yes, will this also change the non-bonded parameters for water - water interaction? 2) Is there a way to add a new ifdef perhaps such that a modified sigma and epsilon can be used for water-protein interactions and the unmodified parameters can be used for water-water interactions? Nonbonded interactions are calculated during the simulation by applying the combination rules defined by the force field. There is no simple way to do this with an ifdef, since that is just in the topology. You can't conditionally apply nonbonded parameters. That just sounds like a recipe for breaking a force field. Not quite right. Parameters for VDW are calculated from the combination rules from the atom-specific values given in [atomtypes] only as a last resort. [nonbond_params] are used in preference to such. So modified protein-water VDW interactions can be introduced by defining all relevant protein atom-TIP3P oxygen [nonbond_params] terms. It may be simpler to modify the [atomtypes] to generate the modified VDW from the combination rule, and introduce the normal TIP3P oxygen-oxygen interaction via [nonbond_params]. Mark, could you please elaborate the method? Mark -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- Quaerendo Invenietis-Seek and you shall discover. -- gmx-users
Re: [gmx-users] (no subject)
On 28/10/2010 2:24 AM, Nilesh Dhumal wrote: Thanks a lot. I made some changes in topology file for rerun simulation. Still I am getting same potential energy. If topology files are different then why the potential energy is same? Occam's razor says that the topologies are not (meaningfully) different. Check your command lines, and compare your .tpr files with gmxcheck. Mark Nilesh On Tue, October 26, 2010 4:29 pm, Mark Abraham wrote: - Original Message - From: Nilesh Dhumalndhu...@andrew.cmu.edu Date: Wednesday, October 27, 2010 1:04 Subject: Re: [gmx-users] (no subject) To: Discussion list for GROMACS usersgmx-users@gromacs.org I used the same .tpr file. I added -dlb yes and - reprod yes during mdrun with rerun option. Still I am not geting why energy, temp, pressure are changing since I have same topology file and .trr file. As I said last time, a parallel rerun cannot reproduce a run unless they're both run under the same conditions. Changing the options for the rerun cannot achieve this, because the original run probably had dynamic load balancing. Mark Is there any bug in rerun option? Nilesh On Mon, October 25, 2010 8:50 pm, Mark Abraham wrote: - Original Message - From: Nilesh Dhumalndhu...@andrew.cmu.edu Date: Tuesday, October 26, 2010 10:50 Subject: Re: [gmx-users] (no subject) To: Discussion list for GROMACS usersgmx-users@gromacs.org I run a test simulation for -rerun. I didn't change the topology file. grompp -f md.mdp -c solvent-bmi-pf6-128.pdb - p solvent-bmi-pf6-128.top -o 3.tpr mpirun -machinefile cp -np 8 mdrun -s 3.tpr -o 3.trr -c solvent-bmi-pf6-128.pdb -e 3.edr -g 3.log with -rerun grompp -f md.mdp -c solvent-bmi-pf6- 128.pdb - p solvent-bmi-pf6-128.top -o 6.tpr mpirun -machinefile cp -np 8 mdrun -s 6.tpr -o 6.trr -rerun 3.trr -c solvent-bmi-pf6-128.pdb -e 6.edr -g 6.log I calculate the total energy by g_energy -f 3.edr -o 3.xvg g_energy -f 6.edr -o 6.xvg The total energy varies between +- 30.00 KJ/mol. It should be constant since I using same topology file and trajectory. Why the total energy is not constant. Those .tpr should be identical - but you can check that with gmxcheck. Reruns do neighbour-searching every step, whereas your normal simulation followed the nstlist setting. That's part of why my earlier advice suggested doing reruns for each simulation you wish to compare. You should be able to get good/better agreement for steps where nstlist directed neighbour-searching in the original run. Also, whether or not constraints have been applied (and when!) could influence the energies to about this degree. I don't recall the details here. Even once you've removed all algorithm-specific sources of difference, there are other sources of non-reproducibility, such as the assignment of particles to DD cells. Your original mdrun probably used dynamic load-balancing, and that cannot be reproduced in the DD used by the rerun. (Or indeed by a repeat of your original mdrun!) Setting -dlb no in the original simulation might be enough to get agreement here, or maybe mdrun -reprod will be required. Mark NIlesh On Thu, October 21, 2010 10:24 am, Mark Abraham wrote: On 22/10/2010 1:17 AM, Nilesh Dhumal wrote: I am doing solvation dynamics for my system. I have system with diatomic (PA---NE)solute surrounded by water molecules. I want to run simulation with two differcent cases. 1. PA charge=0 and NE charge=0 : No charge on solute 2. PA charge=+1 and NE charge=-1 : Charge on solute I want to calculate the energy at each step keeping the solvent configration same. IF I start a simulation with no charge on solute (case1), I have the energy for 1 step. I want to calculate the energy with charge on solute (case 2) with same configration water molecules. Each step I want to calculate the energy with and without charge on solute since the configration of solvent will be same for that step. I was thinking two make two topologies file with charge and with out charge on solute. I don't know how to use them simultaneously during the simulation. Well, you don't use them simultaneously. You run a simulation on whatever you think will generate a relevant conformational ensemble. Then you want to use mdrun -rerun twice on the resulting trajectory, using .tpr files based on .top files corresponding to the two cases in order to create your comparison. 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/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search
Re: [gmx-users] Problems for checking the equilibrium of Ionic liquid system
-BEGIN PGP SIGNED MESSAGE- Hash: SHA1 On 10/28/2010 06:22 AM, hui sun wrote: Dear GMX users, Now, I am doing some simulations about ionic liquids. The problem I am encountering is how to check the equilibrium of Ionic liquid system. The boss asks me to do the RSMD to check it, but I think it is not right. Could you give me some useful advice? Thank you very much for your attentions! With best wishes! Hui Sun Hello, be careful when simulation ionic liquids. I do not know what kind of force field you use, but I assume the one of Lopes and Padua. This force field has problems in reproducing correct dynamics, in other words the diffusion is too slow. We tested this and also other force fields, but dynamics seems to be always the problem. Static properties are represented quite well, however it takes a look time to get into equilibrium. We performed simulations of 100ns and started the analysis after 20ns, if I remember correctly. A possiblity to check if the system is in equilibrium would be to calculate the diffusion in x,y, and z direction seperately and when all plots show normal diffusion, equilibrium should be reached. For further information check my homepage, where you can find some references for papers dealing with this topic. /Flo - -- Florian Dommert Dipl.-Phys. Institute for Computational Physics University Stuttgart Pfaffenwaldring 27 70569 Stuttgart Phone: +49(0)711/685-6-3613 Fax: +49-(0)711/685-6-3658 EMail: domm...@icp.uni-stuttgart.de Home: http://www.icp.uni-stuttgart.de/~icp/Florian_Dommert -BEGIN PGP SIGNATURE- Version: GnuPG v1.4.10 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/ iEYEARECAAYFAkzJNmAACgkQLpNNBb9GiPnauwCg1qw/ul8GyWCXGYp4a9UtDKgn 7r8AoNCbLZKRhXQ5yIgZAkU8AvqroQ22 =nmIy -END PGP SIGNATURE- -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Fwd: RMSF = still confused ?
Hi Lin,
Fair questions.
Worst case: imagine a straight path from A to B, chopped into n=10 stretches d:
A--B
Something goes from A to B linearly and you make a hundred
observations of the position, ten at each stretch. For each stretch i,
the average position is i*d-d/2 and the fluctuation is calculated
about that position, but only includes deviations between -d/2 and
+d/2. The total fluctuation is around (A+B)/2, and includes deviations
with sizes up to +/- (B-A)/2. No surprise that the total RMSF in this
case is far larger than the per bin RMSF. In this example the RMSF
values would end up similar for each bin. If there are some obstacles
along the way, the motion won't be uniform and there will be
differences between bins, but as long as there's a trend going from A
to B (non-equilibrium situation), the total RMSF will be larger than
the RMSF per bin.
Hope this answers your questions.
Cheers,
Tsjerk
-- Forwarded message --
From: Chih-Ying Lin chihying2...@gmail.com
Date: Thu, Oct 28, 2010 at 8:28 AM
Subject: RMSF = still confused ?
To: Tsjerk Wassenaar tsje...@gmail.com
Hi Tsjerk:
Well i am still confused about the RMSF.
1. Consider a particle, it has been observed appearing at 100 positions.
We have ten observations for one time zone. It took the duration
of 10 time zones to get 100 positions
2. The reference position is at the origin. (0,0,0)
3. For each time zones (10 positions), we calculate its standard deviation.
= Then, I can get 10 different standard deviations for 10 time
zones, SD1, SD2, SD3... SD10
4. Further, consider the 10 time zone as a big-whole time zone.
I calculate the standard deviation for the 100 positions, say SD-whole.
5. Then I found
SD-whole SD1,
SD-whole SD2
SD-whole SD3
SD-whole SD4
SD-whole SD5
SD-whole SD6
SD-whole SD7
SD-whole SD8
SD-whole SD9
SD-whole SD10
6. I suppose that
min(SD1~SD10) = SD-whole = max(SD1~SD10) ?
As you said,
the mean squared deviation of the path from A to B around (A+B)/2,
7.
Finally, I test the md-results downloading from your website.
the full set of results
http://nmr.chem.uu.nl/~tsjerk/course/md-tutorial/
THere are 4000 ps. For every 400 ps, I calculate RMSF.
Also, I calculate RMSF(from 0 to 4000 ps) as whole,
I found the same thing that
For each alpha-carbon,
SD-whole SD1,
SD-whole SD2
SD-whole SD3
SD-whole SD4
SD-whole SD5
SD-whole SD6
SD-whole SD7
SD-whole SD8
SD-whole SD9
SD-whole SD10
WHY ???
THank you
Lin
Hi,
Lin, please think your questions over thoroughly in stead of flushing
every thought right to the mailing list. It also helps to stick to a
certain subject (reply) to make sure everything ends up in the same
thread. Maybe it's not a bad idea to read over
http://www.catb.org/esr/faqs/smart-questions.html (again).
Anyway, this question's quite valid :), and Marks answer is a bit
off... The reference structure for computing the RMSF is only used for
fitting; you'd also have seen this behaviour if you'd taken an average
structure from an equilibrium simulation. The fluctuations are around
the mean structure. The larger fluctuations in the first part of the
simulation are caused by the relaxation from the starting structure.
oversimplificationIf you think of the starting structure and the
'equilibrium structure' as position A and B, then the first part shows
you the mean squared deviation of the path from A to B around (A+B)/2,
whereas the next parts give you the mean squared deviation around
B/oversimplification. Note that this gives you a means to assess
whether you've reached equilibrium, albeit not a sufficient measure.
Hope it helps,
Tsjerk
On Tue, Oct 26, 2010 at 8:24 AM, Mark Abraham mark.abra...@anu.edu.au wrote:
If the reference structure from which the fluctuations are measured is taken
from the -s file (check the documentation or code), then a non-equilibrium
trajectory could well have this property w.r.t. a crystal structure.
Mark
- Original Message -
From: Chih-Ying Lin chihying2...@gmail.com
Date: Tuesday, October 26, 2010 15:55
Subject: [gmx-users] g_rmsf = average over # of time frames ???
To: gmx-users@gromacs.org
Hi
From source code = gmx_rmsf.c
g_rmsf computes the root mean square fluctuation (RMSF, i.e. standard ,
deviation) of atomic positions ,
if (devfn) {
/* Calculate RMS Deviation */
for(i=0;(iisize);i++) {
aid = index[i];
for(d=0;(dDIM);d++) {
rmsd_x[i][d] += sqr(x[aid][d]-xref[aid][d]);
}
}
}
count += 1.0;
rmsf[i] = (rmsd_x[i][XX]+rmsd_x[i][YY]+rmsd_x[i][ZZ])/count;
Therefore, g_rmsf is the average of structure deviation over the time
frames.
However, I issued the commands ( C-alpha is selected )
g_rmsf -f abc.xtc -b 0 -e 100 -s abc-crystal.tpr -o RMSF-abc-0th-100th.xvg
g_rmsf -f abc.xtc -b 100 -e 200 -s
[gmx-users] energy group exlusions for frozen structures-use of maxwarn -1?
Hello. I am using a porous structure (MOF). I am interested in how guest molecules inside the porous structure interact with each other and the pores of the structure. As my structure represents a rigid crystalline framework it is usually kept rigid/frozen and only non bonded parameters (LJ and electrostatics) are taken into account for this structure. (I think) I read somewhere on the forum the energy group exclusions should be applied for frozen atoms (also if I don?t do this I get some v. large energies). So I use the following in my .mdp file ; Selection of energy groups energygrps = MOF ; ENERGY GROUP EXCLUSIONS ; Pairs of energy groups for which all non-bonded interactions are excluded energygrp_excl = MOF MOF However on running grompp I get the following warning: WARNING 1 [file MCM_TDI.top, line 32874]: Can not exclude the lattice Coulomb energy between energy groups Why is this? My frozen structure has partial charges but is overall charge neutral. To get around this I have to use maxwarn -1. I am a bit wary of using maxwarn. Can someone advise whether it is sensible to use maxwarn in this case? 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/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Dihedral energy map along two eigenvectors to actual conformations
Found the answer, thanks to personal communication by Dr. Yuguang Mu. Here it is for anyone interested. Each point on the energy landscape corresponds to a time. This time can be tracked simply by using g_anaeig and using -proj to assign a file name and -first, -last to assign eigenvectors of interest. This will generate a file that will show the values of the eigenvectors as a function of time. Subsequently, one can find the time at which the eigenvector values occurred by finding them as a point on the landscape and matching it to the eigenvectors vs time file. Hope my convoluted thinking made sense. :) Cheers, Ali -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] who has the force field file of gromos45a4
Dear all, Right now I plan to use the Gromos 45a4 force field to perform a simulation of carbohydrates. However, this force field is absent from the version of Gromacs (both Gomacs v4.0.7 and v 4.5.1) I'm using. If some one has this force field files, could you, please, pass to me? Best, Liang-- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] who has the force field file of gromos45a4
I answered your same question a few days ago on the mailing list. Please search for my answer if you missed it. Cheers Tom Hong, Liang wrote: Dear all, Right now I plan to use the Gromos 45a4 force field to perform a simulation of carbohydrates. However, this force field is absent from the version of Gromacs (both Gomacs v4.0.7 and v 4.5.1) I'm using. If some one has this force field files, could you, please, pass to me? Best, Liang-- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- Dr Thomas Piggot University of Southampton, UK. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Modifying gromacs
Hi, I want to modify the contribution to Interaction energy of different groups - say I have groups A and B and I want the energy to be scaled as E_AA + 0.1E_AB + 0.5*E_BB. Interaction parameters of each of these groups are set by a forcefield Now after multiple correspondence on the gromacs list I have concluded that there are 3 ways of doing this: 1. Using tables - for this I would have to list non-bonded parameters for all atoms such that the combination rule and the table-potential is used. For the table for BB interaction, scale the Coloumb and VdW interactions in the tables by a factor of 0.5 and so on... However, since tables would have to be supplied for pairs too (tablep), it may not be accurate to supply 6-12 tables with coulomb potential for these pairs. I am using CHARMM and the 1998 paper on CHARMM says that in some specific cases the 1-4 interactions many be scaled which makes me doubt this approach. 2. Forcefield parameters - By defining scaled [nonbonded_params] for all relevant atoms. This will change the VdW interactions, but not sure about the Coulomb interactions. 3. Modifying gromacs - by passing a parameter lambda to gromacs which scales the force/potential by a factor lambda when gromacs calculates force/potential. For implementing option 3, which programs in the gromacs package would be the bes tstarting points for editing the energy contributions of different groups/atoms? As a more general question, how does one run a generalized Hamiltonian REM on gromacs? Pooja -- Quaerendo Invenietis-Seek and you shall discover. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Modifying gromacs
On 29/10/2010 10:14 AM, Sai Pooja wrote: Hi, I want to modify the contribution to Interaction energy of different groups - say I have groups A and B and I want the energy to be scaled as E_AA + 0.1E_AB + 0.5*E_BB. Interaction parameters of each of these groups are set by a forcefield Now after multiple correspondence on the gromacs list I have concluded that there are 3 ways of doing this: 1. Using tables - for this I would have to list non-bonded parameters for all atoms such that the combination rule and the table-potential is used. For the table for BB interaction, scale the Coloumb and VdW interactions in the tables by a factor of 0.5 and so on... Sure. I think that for the above example, you'd need only 2 (maybe 3) table files. However, since tables would have to be supplied for pairs too (tablep), it may not be accurate to supply 6-12 tables with coulomb potential for these pairs. I am using CHARMM and the 1998 paper on CHARMM says that in some specific cases the 1-4 interactions many be scaled which makes me doubt this approach. Yeah, something extra will be required here. Obviously, test and develop on a small toy system that you can compute by hand in a spreadsheet. 2. Forcefield parameters - By defining scaled [nonbonded_params] for all relevant atoms. This will change the VdW interactions, but not sure about the Coulomb interactions. The Coulomb interactions are based on atomic charges, and there's no ready way to scale them differently for different interactions. 3. Modifying gromacs - by passing a parameter lambda to gromacs which scales the force/potential by a factor lambda when gromacs calculates force/potential. For implementing option 3, which programs in the gromacs package would be the bes tstarting points for editing the energy contributions of different groups/atoms? I can think of two ways of approaching this. Efficiency requires that you use the energy group mechanism for your respective groups. Then you need to identify the generated non-bonded lists and do the necessary book-keeping to allocate scale factors to them. Then, either a) pass the scale factor all the way into the non-bonded kernels and use them suitably there, or b) create a copy of the force and energy arrays for each required scale factor, pass matching arrays into the appropriate kernels, and do the scaling on the respective arrays after all kernel contributions have been made, and then add the corresponding array elements. Either way, you'll need a sound understanding of the code in src/gmxlib/nonbonded.c (and the kernels it calls, and how its data structures get created in the neighbour-searching). Make a test system, like a dipeptide in a small box of water, with energy groups, and use a debugger to see how things work. a) is the easiest to develop. The kernels already have a user data pointer you can use, and it would be a fairly simple matter to adapt the generic C kernels to do your scaling. You will give up a lot of performance, however, unless you are prepared to adapt the hardware-optimized kernels similarly (not advised). b) will use somewhat more memory, have a smaller code-change footprint, keep essentially the same performance As a more general question, how does one run a generalized Hamiltonian REM on gromacs? You can't. GROMACS REMD requires a normal MD Hamiltonian supplied in a .tpr. 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/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Modifying gromacs
Thanks Marks! That was very informative. When I try any/all of these I would post the results. Pooja On Thu, Oct 28, 2010 at 10:57 PM, Mark Abraham mark.abra...@anu.edu.auwrote: On 29/10/2010 10:14 AM, Sai Pooja wrote: Hi, I want to modify the contribution to Interaction energy of different groups - say I have groups A and B and I want the energy to be scaled as E_AA + 0.1E_AB + 0.5*E_BB. Interaction parameters of each of these groups are set by a forcefield Now after multiple correspondence on the gromacs list I have concluded that there are 3 ways of doing this: 1. Using tables - for this I would have to list non-bonded parameters for all atoms such that the combination rule and the table-potential is used. For the table for BB interaction, scale the Coloumb and VdW interactions in the tables by a factor of 0.5 and so on... Sure. I think that for the above example, you'd need only 2 (maybe 3) table files. However, since tables would have to be supplied for pairs too (tablep), it may not be accurate to supply 6-12 tables with coulomb potential for these pairs. I am using CHARMM and the 1998 paper on CHARMM says that in some specific cases the 1-4 interactions many be scaled which makes me doubt this approach. Yeah, something extra will be required here. Obviously, test and develop on a small toy system that you can compute by hand in a spreadsheet. 2. Forcefield parameters - By defining scaled [nonbonded_params] for all relevant atoms. This will change the VdW interactions, but not sure about the Coulomb interactions. The Coulomb interactions are based on atomic charges, and there's no ready way to scale them differently for different interactions. 3. Modifying gromacs - by passing a parameter lambda to gromacs which scales the force/potential by a factor lambda when gromacs calculates force/potential. For implementing option 3, which programs in the gromacs package would be the bes tstarting points for editing the energy contributions of different groups/atoms? I can think of two ways of approaching this. Efficiency requires that you use the energy group mechanism for your respective groups. Then you need to identify the generated non-bonded lists and do the necessary book-keeping to allocate scale factors to them. Then, either a) pass the scale factor all the way into the non-bonded kernels and use them suitably there, or b) create a copy of the force and energy arrays for each required scale factor, pass matching arrays into the appropriate kernels, and do the scaling on the respective arrays after all kernel contributions have been made, and then add the corresponding array elements. Either way, you'll need a sound understanding of the code in src/gmxlib/nonbonded.c (and the kernels it calls, and how its data structures get created in the neighbour-searching). Make a test system, like a dipeptide in a small box of water, with energy groups, and use a debugger to see how things work. a) is the easiest to develop. The kernels already have a user data pointer you can use, and it would be a fairly simple matter to adapt the generic C kernels to do your scaling. You will give up a lot of performance, however, unless you are prepared to adapt the hardware-optimized kernels similarly (not advised). b) will use somewhat more memory, have a smaller code-change footprint, keep essentially the same performance As a more general question, how does one run a generalized Hamiltonian REM on gromacs? You can't. GROMACS REMD requires a normal MD Hamiltonian supplied in a .tpr. 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/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- Quaerendo Invenietis-Seek and you shall discover. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
