Thank you very much for the quick reply. Without using tpbconv and mdrun -rerun, I'm not sure how to get the solute-solute energies from gromacs (please let me know it there's a way to do this) -- since there are many more solvent molecules, the total potential energy is dominated by solvent-solvent and solvent-solute interactions. However, the total potential energy from a simulation (not a free energy simulation) with the solute charges set to zero is different than that obtained from a free energy simulation with lambda=1 -- the differences are fairly small, but big enough to be worrisome.
I'll run a system with fewer solvent molecules to get a better handle on what's going on -- thanks again for your suggestions. It seems most people do free energy calculations with a single solute molecule -- in your opinion, is it very unlikely that there is any problem in the way solute-solute interactions are handled in the free energy simulations? Thanks again, Jeff On Wed, Feb 19, 2014 at 7:34 PM, Justin Lemkul <jalem...@vt.edu> wrote: > > > On 2/19/14, 1:05 PM, Jeff Chen wrote: > >> Dear All, >> >> I am trying to calculate the change in free energy that occurs when >> charging multiple type B molecules in an A-B mixture. For example, I have >> a >> mixture of 900 molecules of type A and 100 of type B. When coulomb >> interactions are on, B is identical to A (something like TIP4P water), but >> when coulomb interactions are off, each B is simply a Lennard-Jones >> particle (in the topology file, I've labelled B differently and called the >> molecule SOLB, but all interaction information is identical to A). To do >> the free energy calculations, I perform different simulations at each >> lambda point. For example, for a simulation with lambda=0, I have the >> following in my .mdp file: >> >> free_energy = yes >> init_lambda = 0.0 >> delta_lambda = 0 >> foreign_lambda= 0.05 >> couple-lambda0=vdw-q >> couple-lambda1=vdw >> couple-moltype=SOLB >> >> and when lamba=1, >> >> free_energy = yes >> init_lambda = 1.0 >> delta_lambda = 0 >> foreign_lambda= 0.95 >> couple-lambda0=vdw-q >> couple-lambda1=vdw >> couple-moltype=SOLB >> >> I expected that when lambda=0, B molecules would behave identically to A >> molecules with the coulombic energy fully turned on >> (couple-lambda0=vdw-q), >> and when lambda=1, B molecules would behave like LJ particles, with no >> coulombic interactions, either with each other, or with A molecules. >> >> However, when I look at the coulombic energy for B-B interactions (I used >> tpbconv -o to make a new *.tpr for only B molecules, then used mdrun >> -rerun >> with a trajectory containing only information about B molecules), I find >> that both the short-range and reciprocal coulombic energies are NOT zero >> as >> I had expected for lambda = 1, and are in fact exactly equal to the >> negative of the respective coulombic energies obtained using the same >> trajectory data when lambda=0. >> > > What are the energies without using tpbconv and mdrun -rerun? In other > words, it's necessary to figure out whether this is a problem in the > original simulation with the free energy code (unlikely, given how simple > the system is) or a problem using mdrun -rerun in conjunction with the free > energy code (far more likely). > > -Justin > > -- > ================================================== > > Justin A. Lemkul, Ph.D. > Ruth L. Kirschstein NRSA Postdoctoral Fellow > > Department of Pharmaceutical Sciences > School of Pharmacy > Health Sciences Facility II, Room 601 > 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.
