Re: namdEnergy gives significantly different values from those in log file?

From: Steven Samuel Plotkin (steve_at_physics.ubc.ca)
Date: Mon Nov 24 2008 - 01:50:02 CST

Hi Peter,

I think you are right, and this was why I took the total energies
initially, which includes both KE and PE of the whole peptide plus
water. This is what governs the original sampling probabilities, so I
have discarded the idea of cherry-picking parts of the system. To find
the re-weighted probabilities for each state I have been simply
reweighting each state by the factor:
exp( E_i ( beta_i - beta_targ) )
where E_i is the total energy of the state, beta_i is 1 over kB times
the temperature at that sample time, and beta_targ is 1 over kB times
the target temperature. This should be fine I think, at least for const
volume simulations.

After doing this re-weighting however, the new exponents in the
Boltzmann factors are (when sorted from highest)
       1673.3
         1657
         1638
       1623.2
       1616.6
       1606.1
       1596.6
        .
        .
        .
so one particular state ends up having all the Boltzmann weight in the
partition function (i.e. exp(16.3) more weight than the next weightiest
state). I didn't think this was reasonable, it says entropy isnt
important in the partition function. What are the solutions for this
problem? Perhaps more sampled configurations, but I have 8000 now which
I thought should be enough to at least get a rough idea what the
potential looks like. Maybe I need a narrower temperature range- the
beta_i - beta_targ factors are making the exponents too large. However
previous studies of REMD seem to have gone to 500K without difficulty...
I must be missing something.

I did not add a pressure correction to the replica exchange scripts. I
hadn't thought this through carefully enough and didn't know the REMD
scripts had not accounted for it. According to the reference you
mentioned, there is an extra term in the transition matrix, and an extra
term in the probabilities ~ Vi^N exp(-beta P V_i ). I will have to check
how large this factor is. I don't think the volume fluctuations should
be very large but I could be wrong. Then I would have to rerun my
simulations at constant volume. I assume from reading some previous
posts on the forum that removing all statements pertaining to pressure
and variable volume such as

useGroupPressure yes ;# needed for rigidBonds
useFlexibleCell no
useConstantArea no
langevinPiston on
langevinPistonTarget 1.01325 ;# in bar -> 1 atm
langevinPistonPeriod 100.
langevinPistonDecay 50.

from the conf file should do the trick.

Thanks for any help,
Steve

Peter Freddolino wrote:
> By the way, regarding the text below... I do not think that you should
> be taking the energies of a subset of your system to use in the
> analysis. You really ought to be using the total energy of your system
> for doing the type of analysis that you're doing; otherwise you're
> basically cherry-picking part of the total energy of your system
> without regard for the system as a whole. There are several variants
> of WHAM for analyzing replica exchange data in the literature that
> should be applicable to your problem. Also, the namd replica exchange
> scripts don't contain a pressure correction (Chem. Phys. Lett.
> 335:435439.), so you should be running each simulation at constant
> volume unless you've added it...
> Best,
> Peter
>
> Steven Samuel Plotkin wrote:
>>
>> Dear NAMD users,
>>
>> I have run namdEnergy on a dcd file, printing out All energies. First
>> I noticed the kinetic energy was not printed, then I compared the
>> "TOTAL" potential energy from namdEnergy with the TOTAL energy minus
>> the KINETIC energy in my log file. The two are quite different --
>> about 370 kcal/mol on average. I'm confused by this, I'm not sure
>> what is being calculated differently? I've provided links to my dcd,
>> log, namdEnergy settings, and conf file are below. I've also posted a
>> long version of my question below to help put it in context.
>>
>> (I'm using namd 2.6 and vmd 1.8.6).
>>
>> Thanks for any help,
>>
>> Steve P.
>>
>>
>> FILES:
>> log and dcd file:
>> http://www.physics.ubc.ca/~steve/NAMD/Repx/fold_dse2.job0.0.log
>> http://www.physics.ubc.ca/~steve/NAMD/Repx/fold_dse2.job0.0.dcd
>>
>> Conf files:
>> http://www.physics.ubc.ca/~steve/NAMD/Repx/fold_dse2.conf
>> http://www.physics.ubc.ca/~steve/NAMD/Repx/dse2_base.namd
>>
>> Screenshot of my namdEnergy gui, with settings:
>> http://www.physics.ubc.ca/~steve/NAMD/Repx/Screenshot-NAMDEnergy.png
>>
>> xsc file (if needed):
>> http://www.physics.ubc.ca/~steve/NAMD/Repx/fold_dse2.job0.0.xsc
>>
>>
>>
>> /Long version of my question:/
>>
>> I have been able to modify the Alanin replica exchange tcl code to
>> run on a small peptide which has been solvated and is in ionic
>> solution. I am using namd 2.6 on an x86_64 amd cluster. I ran 16
>> replicas on 8 CPUs across 2 nodes of a cluster, with temperatures
>> ranging from 280K to 440K. My conf file has some new information not
>> in the alanin example because of my solvated box: I have
>> cellBasisVectors, wrapWater, wrapAll, PME, and Pressure control
>> options including langevinPiston. I was careful to choose the
>> langevinPistonTemp to be $NEWTEMP in replica_exchange.tcl.
>>
>> I want to construct a free energy profile for the small peptide.
>> However this must be obtained at one (target) temperature, so I
>> reweighted all states by what their Boltzmann weights should have
>> been at the target temperature. This doesn't work for me: because the
>> system includes the peptide plus the whole water box, after
>> re-weighting, one particular state ends up having all the Boltzmann
>> weight in the partition function. I didn't think this was feasible.
>>
>> So then what I must use is the total energy (KE +PE) of the protein
>> itself, which includes self interactions /plus/ interactions with
>> water and ions. Then I can reweight these energies at the target
>> temperature and get a new partition function and free energy profile.
>> I thought namdEnergy would be the tool for this.
>>
>> As a sanity check, I ran namdEnergy on one of my dcd files, printing
>> out 'All' energies. First I noticed the kinetic energy was not
>> printed, then I compared the "TOTAL" potential energy with the TOTAL
>> energy minus the KINETIC energy in my log file. The two are quite
>> different -- about 370 kcal/mol on average. However the average total
>> potential energy is about 4460 kcal/mol, so the mean error is about
>> 8%. Am I doing something wrong here? Or is this to be expected?
>>
>> At the end of the day what I believe I want is simply the namdEnergy
>> numbers for selection 1= PROTEIN , selection 2= "everything else",
>> added to the PROTEIN-PROTEIN namdEnergies. However the errors present
>> in my sanity check are giving me pause that the namdEnergies are
>> trustable.
>>
>>
>
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:50:08 CST