Hi,
Remember, the 'energy' command is giving you the _total_ energy of
what you provided in the mask, while the 'lie' command is giving you
the _interaction_ energy between the atoms in the first mask and the
atoms in the second mask. If you want to calculate the interaction
energy using the 'energy' command you actually need 3 calculations:
1) The total system energy (energy of components 1 + 2), E12
2) The energy of component 1 (i.e. self energy of 1), E1
3) The energy of component 2( i.e. self energy of 2), E2
Then the interaction energy is the total system energy minus the self
components:
Einteraction = E12 - E1 - E2
I tested that this works with a short test case. I have a FtuFabI
system with NAD and triclosan. I strip off the triclosan and calculate
the interaction energy between FtuFabI and NAD, first with LIE, then
manually with 'energy':
parm FtuFabI.NAD.TCL.parm7
loadcrd FtuFabI.NAD.TCL.nc 1 1 name MyCrd
# Remove triclosan (third molecule)
crdaction MyCrd strip ^3
# Interaction energy via lie
crdaction MyCrd lie ^2 ^1 out temp.lie.dat cutvdw 999 cutelec 999
# Interaction energy via energy
crdaction MyCrd energy E1 ^1 out temp.e1.dat
crdaction MyCrd energy E2 ^2 out temp.e2.dat
crdaction MyCrd energy E12 ^1-2 out temp.e12.dat
Interaction = E12[total] - E1[total] - E2[total]
writedata Interaction.dat Interaction
(The files are located in the $CPPTRAJHOME/test directory if you want
to try it yourself).
Results seem ok accounting for some roundoff etc:
[droe.de6 test] (master)$ cat Interaction.dat
#Frame Interaction
1 -146.8447
[droe.de6 test] (master)$ cat temp.lie.dat
#Frame LIE_00002[EELEC] LIE_00002[EVDW]
1 -91.6226 -55.3213
[droe.de6 test] (master)$ echo "-91.6226 + -55.3213" | bc -l
-146.9439
Let me know if that clears things up.
-Dan
On Tue, Jun 24, 2025 at 9:59 AM Игорь Упоров <iuporov.enzyme.chem.msu.ru> wrote:
>
> Daniel,
>
> Thank you very much for your quick reply.
>
>
>>
>> The manual is your friend here. The 'lie' command by default uses a
>> cutoff for non-bonded interactions:
>>
>> "Electrostatic and van der Waals interactions will be calculated for
>> all atom pairs. A separate electrostatic and van der Waals cuto can
>> be applied, the default is 12.0 Angstroms for both."
>
>
> It's written in an ambiguous way. I've understood this text like - on default calculations for all atom pairs, not like cutoff radii 12.0 Angstroms on default.
>>
>>
>> The "simple" non-bonded energy model used by the 'energy' command by
>> default has no cutoff:
>>
>> "Note that the non-bonded energy terms for 'simple' do not take into
>> account periodicity and there is no distance cut-o ff. "
>>
>> I suspect this is where the discrepancy lies.
>>
> OK. I've calculated using
>
> lie :330 :1-329 out Enbond1.csv nopbc cutvdw 30 cutelec 30
>
> Diameter of my complex is 60A, so a cutoff 30A is enough.
> Results are:
>
> #Frame LIE_00001[EELEC] LIE_00001[EVDW]
> 1 -81.7649 -31.6603
>
> Better than used to be, but far away from results of "energy" command
>
> Eelec EVDW Total
> -151.7 -35.7 -187.4
>
> VDW terms [A/(30)^12 - B/(30)^6] are extremely small, but we see a difference 4 kcal/mol in VdWaals energy.
> Which command should I use to calculate real interaction energy?
>
> Thank you one more time.
> UIV.
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Jun 24 2025 - 08:30:02 PDT
